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Abstract 



In this thesis we investigate the properties of ultracold Bose gases in one and two dimen- 
sions. Low-dimensional systems have several striking differences from their three-dimensional 
counterparts, most notably the absence of Bose-Einstein condensation (BEC) in the infinite 
homogeneous case. Recent experimental progress has brought low-dimensional systems within 
reach in the laboratory, and we provide numerical simulations with experimentally relevant 
parameters. We present simulations of the Berezinskii-Kosterlitz-Thouless (BKT) phase tran- 
sition in two dimensions and investigate equations of motion for a Bose gas constrained to one 
dimension. 

Recent experiments have attributed superfluidity in two-dimensional systems to the BKT 
phase transition. We perform classical field simulations using the projected Gross-Pitaevskii 
equation (PGPE) formalism to model the two-dimensional Bose gas at finite temperature. We 
confirm the presence of the BKT phase via the observation of two unique features: algebraic 
decay of the first-order correlation function; and vortex pair unbinding. Unbinding of vortex 
pairs at the BKT transition is clearly demonstrated via a coarse-graining procedure which re- 
veals unpaired vortices. The BKT transition temperature identified via correlations and vortex 
behaviour agrees well with the temperature deduced from a calculation of the superfiuid frac- 
tion. Surprisingly, we observed no separation between the temperature of the BEC transition 
— which is present due to the finite simulation size — and the BKT transition. We relate our 
results to experimental observations and show that an interpretation based on BKT physics is 
justified. 

In investigating the two-dimensional system we found it necessary to compute the super- 
fluid fraction. Calculating the superfiuid fraction is a delicate procedure because it involves 
connecting the macroscopic quantum phenomenon of superfluidity to a detailed microscopic 
simulation. We present an efficient method which overcomes this difficulty using a tensor de- 
composition of the momentum density autocorrelations. Our method gives results consistent 
with the other physical properties of the BKT phase in two-dimensional systems. 

One-dimensional effective equations allow for efficient simulation of very elongated sys- 
tems of ultracold Bose gases. We use a Gaussian ansatz and Lagrangian approach to derive 
an effective equation for a Bose gas constrained to one dimension at zero temperature. In 
some respects our method outperforms alternative one-dimensional equations such as the non- 
polynomial Schrodinger equation. However, our scheme is found to be inherently unstable in 
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a wide range of cases. We analyse the instability and find that it is an inherent feature of a 
whole class of nonlinear ansatze. 
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Chapter 1 

Introduction 



In this thesis we investigate the behaviour of dilute, degenerate Bose gases confined to one 
and two dimensions. This introductory chapter is mostly dedicated to explaining the physics 
contained in the previous sentence: The quantum statistics of bosons, the nature of degeneracy, 
and the physical differences between three and lower dimensions. 

It is a surprising fact that elementary particles of the same type are not distinguishable, 
even in principle. The first hints of the importance of this fact appeared even before the advent 
of quantum mechanics as a resolution to the Gibbs paradox of classical statistical mechan- 
ics: The calculated entropy of an ideal gas is not extensive^ unless the atoms are treated as 
indistinguishable. As it turns out, there are exactly two consistent statistical behaviours for 
indistinguishable quantum particles in three dimensions^. Bose discovered the first of these in 
1924 [17] in a successful attempt to derive Planck's blackbody radiation formula from first prin- 
ciples. Bose had invented a new way to count multi-photon states, which Einstein immediately 
applied to atoms to predict the equation of state [36], and in 1925 the "condensation" [37] of 
what became known as the ideal Bose gas. Bose-Einstein condensation (BEG) was invoked in 
the following decades to explain various condensed matter phenomena including superfiuidity 
in liquid helium and low temperature superconductivity. In 1926 Fermi [40] and Dirac [33] 
discovered the second type of quantum statistics, motivated partly by the behaviour of elec- 
trons as newly articulated in the Pauli exclusion principle. Particles obeying Bose-Einstein and 
Fermi-Dirac statistics are now known as bosons and fermions respectively. 

We refiect briefiy on what exactly we mean by "particle statistics". Statistical mechanics 
starts from a probability distribution over the system microstates^ and uses this distribution to 
calculate expected properties of the system. The distribution is specified by some macroscopic 

measured quantity is extensive if it is directly proportional to the size of the system; for a new system 
composed of two exact copies of the system joined together the extensive quantity is doubled. 

^In two dimensions there is in fact a continuum of possibilities known as anyonic statistics. The systems 
considered in this thesis are embedded in real three-dimensional space however, so anyonic statistics are irrelevant 
at the level of the constituent particles. 

■^A microstate is a full microscopic description of the details of a physical system; a state of maximal 
information. 
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variables (for example the temperature) but before the distribution can be given, the possible 
microstates must be described. Bosons, fermions and classical particles differ in exactly which 
microstates are possible. For a noninteracting system such as the ideal gas, the state of the 
full system may be given by specifying which single-particle states, or modes'^, are occupied. 
Here is where the difference arises — any number of bosons may occupy a given mode, while 
the occupation is constrained to zero or one particle for fermions. Labelling the states of the 
full system using only the occupation of modes incorporates particle indistinguishability. For 
classical distinguishable particles, we would also need to keep track of which particles are in 
each mode, in addition to how many. 

Quantum statistical behaviour is connected to intrinsic angular momentum by the celebrated 
spin statistics theorem: particles of half-integer spin (in multiples of h) are fermions while 
particles with integer spin are bosons. Although the elementary constituents of normal matter 
are all fermions, bosonic statistics may be observed in composite particles — atoms for example 
— composed of even numbers of fermions. For this to be true, it is sufficient that the interactions 
between composites be at much lower energy than the internal energy levels [35]. 

In order to observe the consequences of quantum statistics a system must be degenerate. In 
this context "degenerate" means that the number of energetically accessible modes is similar 
to or smaller than the total number of particles; it is only with this kind of crowding of the 
available states that the differences in counting microstates become apparent. Degeneracy may 
be reached by either increasing the density, which increases the number of particles per state, 
or reducing the temperature which reduces the number of accessible states. For the particular 
case of an ideal gas, it is convenient to discuss degeneracy in terms of the de Broglie wavelength: 
At a given temperature T, the expected energy of a classical ideal gas particle is |/c_bT for each 
degree of freedom, according to the equipartition theorem, where ks is Boltzmann's constant. 
In three dimensions we therefore have (p^/2m) = |/cbT, and using de Broglie's formula X = h/p 
gives a value X = h/ y/^mkBT for the typical quantum wavelength. It is conventional^ to define 
the thermal de Broglie wavelength as 



With this definition the number of accessible states can be estimated as, Z = V/ X\^ where V is 



^It is easy to confuse the terminology when referring to "single-partiele states" versus "states of the system as 
a whole" . We hope to avoid this problem by following a common practice from optics where the single-particle 
states are referred to as modes. 

^The convention for Ajb is selected to absorb the constants in the expression for the ideal gas partition 
function. 
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the volume of the system^. Quantum degeneracy occurs when the number of particles is greater 
than or approximately equal to the number of accessible states, that is, N > Z. Rearranging, 
we arrive at the degeneracy condition 

n > Ki (1.2) 

for the density n = N/V, which has a nice physical interpretation: Degeneracy occurs when 
the size of the quantum wavepacket is comparable to or larger than the interparticle spacing. 

Apart from Bose's work, the first historical uses of quantum statistics were to compute 
properties of the ideal gas in the quantum degenerate regime. Nevertheless, there was no 
experimentally accessible system with both weak interactions and strong quantum statistical 
behaviour until 1995 when the first degenerate, dilute Bose gases were created in the laboratory 
[4, 25]. This feat was followed shortly afterwards with the creation of a degenerate Fermi gas in 
1999 [32]. There were several good reasons why it took seventy years for experimental systems 
to catch up with the theory. On the one hand, quantum degeneracy occurs in many sorts 
of condensed matter systems, including striking examples such as superfluidity in liquid "^He. 
However, these systems have strong interactions between particles which renders the simplest 
theories useless from a quantitative perspective. On the other hand, gases with weak interac- 
tions are nowhere near degeneracy at easily accessible temperatures, and tend to form liquids 
and solids when cooled. At low densities the formation of liquids or solids from atomic gases is 
suppressed, and cooling to extremely low temperatures — of the order of 100 nK — allows ex- 
periments to reach the degenerate regime. Perfecting the experimental tools necessary to trap 
and cool atomic clouds to such extraordinarily low temperatures was a major undertaking. 

From now on we focus exclusively on bosons which are the topic of this thesis. The typical 
ultracold boson experiment begins with a gas of neutral atoms — ^^Rb for example — in a 
vacuum chamber. In the first stage the atoms are captured in a magneto-optical trap and cooled 
via laser cooling to temperatures of order 100 fiK [98]. Atoms are then loaded into a tighter 
magnetic or optical dipole trap and further cooled by successively removing the most energetic 
atoms — a process known as evaporative cooling. During this process the gas undergoes 
condensation to form a BEG at temperatures of order 1 fiK to 100 nK, with total number of 
atoms ranging from 10^ to 10^ [96, §1.1]. The trapping potential is usually well approximated 
by a parabola, with the spatial extent of the atomic cloud on the order of 10-100 yU.m, depending 
on the trap anisotropy [5]. 

The precision and control available to experiments has continued to improve over the last 
15 years. Experiments can now control the interaction strength via Feshbach resonances, and 



^Hcrc Z is actually the canonical partition function; if the energy scale is adjusted such that the lowest 
energy state has energy zero. Z counts the accessible states. (See, for example, [115, §6.1].) 
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detailed control of the trapping potential is available using combinations of lasers and magnetic 
fields. Ultimately, the experimental accessibility of weakly interacting and precisely controllable 
quantum gases has been a great resource for basic quantum physics: It has allowed rigorous 
comparison and evaluation of first principles approaches to quantum field theory. 

1.1 Bose statistics and condensation 

The phenomenon of BEC is a phase transition, with a normal gaseous phase at higher tem- 
perature, and a fluid with macroscopic quantum behaviour at low temperature. In the low 
temperature phase, a macroscopic number of particles all occupy the same quantum state, in a 
sense to be made more precise below. For the particular case of the noninteracting homogeneous 
gas in three dimensions, the expected occupation of the ground state is 



where Tc is the transition temperature. We discuss the value of Tc and provide a derivation of 
this basic relation in section 2.1.2. 

Unlike most phase transitions, BEC occurs even in systems without interparticle inter- 
actions: it is a purely statistical effect. To understand qualitatively how this happens, we 
consider the effect of Bose statistics on a two mode toy system containing exactly particles. 
For definiteness, let us think of these modes as a pair of "left" and "right" potential wells. 

First let us consider the case where each well has the same energy, and count the number of 
states available to the system as a whole. For classical distinguishable particles, each particle 
may be in the left or right well independently of the rest, so adding a particle to the system 
multiplies the number of available states by two; the total number of states for particles is 
then 2^. On the other hand, for indistinguishable particles a state of the system as a whole is 
fully specified by listing the number of particles in each well. The number of particles in the 
left well, is between and particles, and is always equal to — N^, so there are 
only A^ + 1 states. 

If both modes have the same energy then all states of the system have the same energy 
and are therefore equally likely''. The probability of having all A^ particles in the left well is 
then 1/2^ for the distinguishable case, but 1/(A^ + 1) for indistinguishable particles. From this 
simple example, we can already see how Bose statistics exponentially enhances the probability 
of finding the system with all particles in a single well, simply because of the way states are 

""The assumption that states with equal energy are equally probable is the founding assumption of statistical 
mechanics, see, for example, Ref. [115, §7.1]. 
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Figure 1.1: Bose statistics versus classical distinguishable statistics in a two mode toy system, 
(a) Count of system states with a given number of particles in the left well for the distinguishable 
case (blue) and bosonic case (green). This is proportional to the probability density p{Nl) when 
A = 0. (b) Probability density p{Nl) for finding the system with Nl particles in the left well, 
when A/ksT = 1 [colours as in (a)]. Both (a) and (b) show statistics for a system of iV = 7 
atoms, (c) Expected fraction of particles in the left well as a function of A/Zc^T for the bosonic 
case with three different atom numbers. The distinguishable case is included for reference and 
is the same for any N. 
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counted. One way to visualise this is to plot the number of states with a particular A^^ as in 
Fig. 1.1(a) — even for the very small value of = 7 shown, there are many more distinguishable 
states where the particles are evenly distributed between the wells. 

If we now let the right well have energy A greater than the left, the probability distribution 
p{Nl) is skewed proportionally to q^^r^/i^bT ^ shown in Fig. l.l(b)^. Nevertheless, the 
distinguishable case may still have a peak at nonzero A''^, due to the statistical pressure toward 
equal numbers; this effect becomes much more pronounced for larger, more realistic values of 
N . Taking the expected value of A'^^. from the probability distribution, we can also calculate^ 
a "condensate fraction" — the expected fraction of atoms in the lowest energy state, {Nl)/N. 
Figure 1.1(c) shows the condensate fraction {Nl)/N, as a function of A/ksT, for various total 
number A^. It is clear that for bosons the condensate fraction rapidly converges to 1 as A^ 
increases, regardless of the size of the energy gap between the wells. On the other hand, for 
distinguishable particles the fraction in the left well is independent of A^. It is worth noting 
that in this toy system — or any system with a finite number of modes — there must always 
be a mode with an extensive population, so strictly speaking the existence of condensation is a 
trivial fact in this case. Even so, it is instructive to observe the differences between the bosonic 
and distinguishable cases. 

An infinite homogeneous gas has an infinite number of modes, so the existence of a nonzero 
condensate fraction in the lowest energy mode is a nontrivial fact in this case (see section 
2.1.2). The existence of BEG in a system with infinitely many modes turns out to depend on 
the density of states at low energies. In the three-dimensional homogeneous case, there are 
"sufficiently few" low-energy modes meaning that particles are forced to gather in the lowest 
energy mode and BEG occurs. However, in one and two dimensions the density of low-energy 
modes is larger and BEG is forbidden, at least in the homogeneous case. We emphasise that 
the addition of a trapping potential modifies the density of states, which can lead to BEG even 
in low dimensions. 

1.2 Low-dimensional systems 

The physical behaviour of one- and two-dimensional systems can be qualitatively different from 
a three-dimensional (3D) system containing the same constituent particles. A dramatic example 
is the homogeneous ideal Bose gas which undergoes the BEG phase transition in 3D, but never 
condenses in lower dimensions at any nonzero temperature. This absence of BEG is a specific 

'^The Boltzniann factor g-^aA/kBT g^j-jggg from analysing the system in the canonical ensemble, see, for 
example, Rcf. [115, Ch. 6]. 

^Thc calculation of {Nl)/N is immediate in the canonical ensemble via direct summation over states. 
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case of the more general Mermin-Wagner-Hohenberg theorem [86, 66], which forbids long-range 
order in a whole class of one- and two-dimensional systems. 

The qualitative differences between three and lower dimensions are due to both geometric 
and topological effects. For example, the existence of BEC depends critically on the density of 
states at low energy, which for the homogeneous case scales as the surface area of a (rf— l)-sphere 
in momentum space, where d is the dimension; this is a geometric effect (see section 2.1.2). On 
the other hand, the nature and existence of elementary topological excitations such as vortices 
also depends on the dimensionality [19, Ch. 9] and these excitations have a significant effect on 
the thermodynamics of the system. 

At first sight one might question the relevance of low-dimensional systems, given that the 
world is really three-dimensional. However, many 3D systems can be made to act as if they 
are low-dimensional by tight confinement. When the confinement is sufficiently strong and the 
temperature low enough, the transverse degrees of freedom are effectively "frozen out" and all 
dynamics occurs within the plane or line. More precisely, when interactions are unimportant it is 
clear that freeze out occurs for fiw ^ ksT, where we have compared one quantum of transverse 
excitation energy f\w with the characteristic magnitude of the thermal fluctuations. This 
estimate holds even when interactions are important, as discussed in Ref. [96, Ch. 15]. 

Ref. [54] presents a useful characterisation of low-dimensional regimes in terms of length 
scales: Let us consider a BEC in a cylindrically symmetric system, where the extent of the 
trapped cloud is given by radii i?_L and Rz in the radial and longitudinal directions respectively. 
As discussed in Ref. [54], there are two additional length scales relevant in this problem, the 
healing length^° ^ and the scattering length a^. In many experiments, we have the situation 
-R2 ^ ^ ^ a^, in which case the dynamics are fully three-dimensional. In contrast, the 2D 
and ID cases correspond to the regimes R±_ ^ C, > Rz and Rz ^ C, > R±, where the BEC 
is conventionally called "pancake-shaped" and "cigar-shaped" respectively. In all these cases 
we assume that R±,Rz ^ ag so that the scattering remains three-dimensional, even when the 
dynamics is two- or one-dimensional. This is called the quasi-low-dimensional regime [54]. 

We briefly discuss some unique properties of two-dimensional systems. As noted above, 
condensation does not happen in a homogeneous 2D system. However, the 2D Bose gas sup- 
ports topological defects in the form of vortices, and in the presence of interactions can instead 
undergo a Berezinskii-Kosterlitz-Thouless (BKT) [9, 75, 104] transition to a quasi-coherent su- 
perfiuid state. By quasi-coherent we mean a state with algebraic^^ decay of spatial correlations, 

^"The healing length ^ = l/^/STmoT is the length scale on which a condensate recovers ("heals") from the 
effects of a local perturbation, where n is the density. ^ may be estimated as the length scale on which the kinetic 
term in the in the Gross-Pitaevskii equation (described in section 2.3) is equal to the size of the interaction 
term far from the perturbation. For details see Ref. [96, §6.4]. 

11 "Algebraic" decay is the name conventionally used in this context for a simple power law decay. 
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rather than the traditional off diagonal long-range order of a coherent condensate [124], or the 
short range exponential decay characteristic of a high temperature phase. This is one of several 
unique features of the BKT phase. 

The most striking feature of the BKT transition is the pairing of topological defects in 
the low temperature phase. In the superfluid case this corresponds to pairing of quantised 
vortices of opposite rotation, which begin to unbind as the temperature is increased through 
the transition; the high temperature phase is a plasma of bound and unbound vortices [75]. 

A simple argument shows why vortex unbinding is related to the destruction of superflu- 
idity above the BKT transition. Consider a channel filled with a superfluid, and transport a 
single unpaired vortex from one boundary of the channel to the other. This imparts a single 
quantum of velocity to the fluid within the channel^^, and provides a mechanism for the decay 
of supercurrents. On the other hand, transport of pairs of vortices does not impart any net 
momentum to the channel. We defer further presentation of the aspects of the BKT transition 
relevant to 2D Bose gas experiments until chapter 3. 

In addition to 2D BKT physics, we also consider BECs in ID systems in this thesis. In 
contrast to the 2D case, we make no particular attempt to investigate features unique to one 
dimension, since our focus is on developing a numerical approximation technique. However, we 
do rely on one uniquely ID feature for testing our numerical method: In ID there are stable 
soliton solutions with well understood behaviour; in higher dimensions such solitons would be 
unstable, and eventually decay via the "snake instability" [3]. 

Experiments in cold gases have been able to access the low-dimensional regime since 2001 
when Gorlitz et al. created both ID and 2D systems using magnetic and optical traps respec- 
tively [54]. The BKT transition was first observed in liquid helium thin films [11], but observing 
aspects of BKT physics in an ultracold Bose gas proved difficult until a series of experiments by 
the Dalibard group at ENS in Paris [121, 60, 77]. In the Paris experiments, an optical lattice 
was used to create two parallel pancake-shaped BECs. Using these as mutual phase references 
allowed unpaired vortices to be detected in an interference experiment [121]. Various aspects 
of BKT physics including the algebraic decay of correlations were subsequently observed in 
Refs. [60, 77]. Further experimental work by other groups has followed [117, 23], further adding 
to the evidence of BKT physics in the 2D cold gas system. 



^^It is easiest to see how the single quantum is imparted by considering a toroidal channel where the superflow 
represents angular momentum, and transporting a vortex from the centre to the outside. 
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1.3 Outline of the thesis 

After the current introductory chapter, some background theory is introduced in chapter 2, 
including the two main pieces of formahsm used later in the thesis: The time dependent Gross- 
Pitaevskii equation (to be used in chapter 5) is introduced in section 2.3; the classical field 
methods for finite temperature calculations (in particular the projected Gross-Pitaevskii equa- 
tion, to be used in chapter 3) are introduced in section 2.4. An attempt has been made to keep 
the background self-contained, starting from the level of an advanced undergraduate. For this 
reason, a section on basic Bose statistics is included, along with a discussion of how nonrela- 
tivistic quantum field theory arises as a generalisation of first- quantised many-body quantum 
mechanics. 

The two-dimensional experiments discussed in the previous section present a challenge for 
theory. One particular issue is to understand the relation between the BKT phase and BEG, 
which exists in real experiments due to the harmonic trapping potential. In chapter 3 we make 
use of classical field methods to simulate the finite temperature physics of 2D systems, with a 
view to understanding this question and to understand which aspects of BKT physics may be 
observed in current experiments. 

As part of this work, we developed a method for calculating the superfiuid fraction from 
a classical field simulation. This applicability of this method is more general than the specific 
use to which it is put in chapter 3, so we describe the derivation in isolation as chapter 4. 

We change focus in chapter 5 to examine effective equations for ID systems. Motivated by 
the desire to study dispersive shock waves in high resolution and perform expansions of suddenly 
untrapped atomic clouds, we consider a variational ansatz for the quasi-lD time evolution. We 
implement the method numerically and characterise the behaviour using several test problems. 

The thesis concludes with a summary in chapter 6. Several appendices are included: ap- 
pendix A describes subtleties surrounding differentiation of Lagrangian functionals with respect 
to complex fields. Appendices B and G describe additional details pertaining to chapters 3 and 
5 respectively. Finally, appendix D contains verbatim copies of Refs. [18, 84]. These two papers 
were produced during the author's PhD candidature, but were not related to the main theme 
of the thesis. 



Chapter 2 

Background Theory 



In this chapter we describe some basic theoretical background followed by the specific theoretical 
methods used in the thesis. We start with some elementary remarks about condensation in Bose 
gases and then show how nonrelativistic bosonic quantum field theory arises from single-particle 
quantum mechanics. We next derive the Gross-Pitaevskii equation from the variational point of 
view. The projected Gross-Pitaevskii equation is covered in some detail, including a discussion 
of the motivation and procedure for taking the classical limit. 

2.1 Bose statistics and BEC in noninteracting systems 

Bose-Einstein condensation is a unique phase transition because it exists even in the absence 
of interactions between particles: it is a purely statistical effect. Without interactions, the 
behaviour of a Bose gas can readily be analysed using the methods of elementary statistical 
mechanics. The following outlines the argument, showing that the occupation of the lowest 
energy state is extensive below some transition temperature, which we then compute. 

2.1.1 The Bose distribution 

In this section we derive the Bose distribution from basic statistical mechanical considerations. 
Before we start it is worth noting that the only quantum mechanical ingredient required is 
indistinguishability of particles. When we talk about microstates of the system, we may think 
about either classical states with a particular total energy and number of particles, or quantum 
eigenstates of the Hamiltonian. 

In the grand canonical ensemble, the probability of observing a given microstate s of a 
system is proportional to the Gibbs factor: 




(2.1) 
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Here E{s) is the energy, N{s) is the number of particles in the microstate and /3 = l/ksT is the 
inverse temperature. The Gibbs factor comes about by considering a system that can exchange 
energy and particles with a much larger reservoir, and applying the fundamental assumption 
of statistical mechanics — that all microstates in a closed system are equally likely — to the 
closed combination of system and reservoir (see, for example, [115, §7.1]). The normalisation 
constant for P is the sum of the Gibbs factors for all states, known as the partition function 

2 = ^e-/3m-MA^W]. (2.2) 

with this definition we have P(s) = Z^^e"^^^^^^~^'^^^^^. 

We consider a system with a number of discrete single-particle quantum states or modes, 
{|0), |1), |2), . . .}. Letting the single-particle energies for the modes be and considering the 
case where particles have no interaction energy, the total energy is E{s) = Xli ^«'^«('^)' with 
corresponding total number N{s) = Yli'^^ii^) where rij is the occupation for the zth mode. 
With these expressions, the partition function for a noninteracting system is 

2 = ^g-/3Eife-M)".W (2.3) 

s 

= ^ JJ" g-/3(ei-M)"i(«) (2.4) 
s i 

where the product ranges over all modes. 

To proceed further, we need to make use of Bose statistics to specify over which states 
the sum occurs. Because bosons are indistinguishable, a state is uniquely defined by the 
number of particles in each mode, s = [ni{s) , n2{s) , . . .]. Therefore, after some thought^ Z can 
be rewritten in the simpler form 

oo 



Z = JJ ^ e-^^^^-'^)" (2.5) 

i n=0 

= l[2,, (2.6) 



^To make this transformation slightly less mysterious, consider only two modes. Enumerating the pos- 
sible states as [0, 0], [1, 0], [0, 1], [2, 0], [1, 1], [0, 2], [3, 0], [2, 1], [1, 2], [0, 3], . . . , and abbreviating a = e-'^^^^-^'\ 
b = e~^^'^^~'^\ we see that we're dealing with the expression 

Z = Y^ a»i(«)6"2(s) = i + a + b + a'^ + ab + b^ + + a% + ab^ + + ... 

S 

OO oo 

= (1 + a + + + . . . )(1 + 6 + 52 + + . . . ) = 51 II 



n=0 n=0 
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where we have defined the single-mode partition function Zi = Yl'^=o^~^^'^^~^^"'- We note that 
the simphcity of Eq. (2.5) is what makes the analytical calculation tractable, and is peculiar to 
the grand canonical ensemble. 

An explicit formula for Zi is easily obtained by summing the geometric series: 



oo oo ^ 

Z. = ^e-""-")" = Y: (e-«..-))" = (2.7) 



n=0 n=0 

This puts US in a position to calculate the expected occupation of the jth mode. Because of 
the factorised form of Eq. (2.5), we can show that functions of Uj are normalised using only Zj 
rather than Z, 

^ oo 

(fin,)) = E fM'))P(') = ^ E fir^)e-^^''-'^'' (2-8) 

s n=0 

for any function /. For the particular case of f{n) = n we can evaluate the sum^, leading to 
the famous Bose distribution function 



(n,) = — ^ . (2.9) 



2.1.2 Condensation in the 3D homogeneous Bose gas 

Having derived the Bose distribution, we turn our attention to the statistical mechanics of Bose 
condensation in a 3D gas. The expected total number of atoms in the system is given by the 
sum over modes 

^^ = E("-> = Ep;5:^- (2-10) 

i i 

For a given system the modal energies {e,} are known, so at a particular total number 
and inverse temperature /3 we may in principle use this equation to determine /i. However, 
solving this equation is not possible except in the simplest cases, and we therefore turn to 
approximations. 

The general procedure is to split the atoms into two groups: the expected numbers in the 
ground and excited states, Nq and Apx respectively. The original argument due to Einstein 
is that in the 3D homogeneous gas there is an upper bound on the number of excited state 
atoms Nf.^ at any given temperature. For any system with more atoms than this bound, the 



^These kinds of sums may be reduced to the sum of a geometric series using the common derivative trick; 

d I 



n=0 n^O ^ n^O 



y— In X 



dy 1 ~ ey 



y— In X 
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remaining atoms must be found condensed in the ground state [37]: 

[. . .] something similar happens as when isothermally compressing a vapour beyond 
the volume of saturation. A separation occurs; apart "condenses" , the rest remains 
a "saturated ideal gas" . 

The expected number of excited state atoms is the sum over excited state modal occupations, 

A^ex = V \ • (2-11) 

To estimate N^^ analytically we need a few transformations and approximations, starting by 
rewriting the sum as a sum over the mode energies rather than state indices. Defining the 
discrete density of states 5'dsc(e) to be the number of modes at energy e, we have 

iVe. = ^^dsc(6) ^,(j)_^ . (2.12) 
e 

Next, we recognise that A^o > implies yU < eo, which in turn implies l/(e^'^'^'~^^ — 1) < 
1/ (^e^i^i-'^o) _ gQ that we can replace fx with eg to bound A^cx from above: 

iVex<E^dsc(6) ^,(j^)_^ . (2.13) 
e 

In the particular case of the ideal gas, eo = and we have 

iVex< J]^?dsc(6)^^. (2.14) 
e 

Unfortunately evaluating this sum is not straightforward, so to progress further we transform 
it into an integral via an approximation. We have 

A^ex< / de gie) ^^^J^^ _ (2.15) 

where g{e) is the density of states that will correspond to a smoothed version of 5'dsc(e)- 

The most elementary method for computing an appropriately smoothed g is to first evaluate 
the total number of states G{e) with energy less than e (see, for example, [96]). With the 
simplest approximate evaluation, G{e) turns out to be smooth and we may define a smooth g 
using g{e) = For the 3D gas in a finite periodic box of size L x L x L, the possible modes 
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are indexed by the integer wavenumber q G Z and have energy 

^ q = 2.16 

The number of modes with energy less than e is then approximately the volume of the sphere 
of radius ^J^rnJJe[2W^ containing the set {q G R'^: -E'(q) < e}, that is, 

47r / mL'^e^ 



Differentiating G leads to the density of states 

, , dG 4V27rL3m3/2 

»w = * = — P — '2-'**) 

which we see varies as the square root of e. 

Putting this into Eq. (2.15) and performing the integral^, we finally arrive at 

K.<y C(|)T^/' (2.19) 

where C{'^) = Xlfcli 1/^" is the Riemann zeta function, with (^(|) 2.61. This expression tells 
us that there must be some critical temperature Tc below which is strictly greater than Nc-x, 
and the remaining — N^^ atoms must go into the ground state. At the critical temperature, 
A^ = Ncx{Tc) which implies 

T, = ^. (2.20) 

The fraction of atoms in the ground state as a function of temperature is 

N - N / 7^ \ 3/2 

Determining the condensate fraction using the standard derivation presented above depends 
critically on the functional form of the density of states. In turn, the density of states depends on 
the dimensionality and any trapping potential, both of which modify the possible single-particle 



■^The integral may be evaluated by expanding in powers of e ^. followed by a change of variables and 
recognising the Riemann zeta function ^ and gamma function F: 

/ dx^— = y^ dxxi/2e-™ = Vn~3/M (ixxi/2e-^=C(i)r(f). 
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quantum states. Indeed, in the two-dimensional homogeneous gas the integral in Eq. (2.15) 
fails to converge, suggesting — though not by itself proving — an absence of condensation. 

2.2 Field theory and interacting Bose gases 

In this section we develop some basics of quantum field theory in a form that applies to cold 
gases of atomic bosons. In particular, the theory developed below is strictly nonrelativistic, 
and assumes bosonic symmetry for the many-particle quantum state. 

2.2.1 Second quantisation 

The formalism of second quantisation^ is an elegant and compact way to express the quantum 
mechanics of identical particles. Second quantisation succeeds by incorporating the symmetries 
of the many-particle state into the theory from the start. At heart, it may be viewed as 
a convenient reformulation of basic many-particle quantum mechanics, extended to allow for 
superpositions of states with different particle number. There are many good explanations of 
second quantisation in the literature; for additional exposition we refer the reader to Refs. [2, 
122]. 

Many-particle symmetries and the number state basis 

To describe second quantisation we start with the quantum mechanics of an particle system. 
In the general case, the Hilbert space for a system of distinguishable particles is the tensor 
product of A^ one-particle Hilbert spaces (see, for example, [93, §2.2.8]): 

'H = Hi®n2®---®'HN. (2.22) 

Accordingly, the possible states of A^ identical particles live in the Hilbert space 

'H = 'Hi®Hi®---®Hi ='H^ . (2.23) 

N times 

States from the natural tensor product basis for l-L^ have the form^ |ji)|j2) ■ ■ ■ \ jN)i where the 
ii are positive integers and we use {|1), |2), . . .} as a basis for "Hi. 

^The name "second quantisation" is somewhat obscure. Historically, one method of deriving a non-interacting 
quantum field theory was to apply the procedure of canonical quantisation to the Schrodinger field (see [122, 
page 81] for a discussion). Insofar as one step of quantisation has already been used to obtain the Schrodinger 
equation, it appeared that this procedure constituted quantising a second time. 

^We omit the tensor product symbols between kets for brevity; \a)\h) = \a) ® \b). 
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It is possible to solve problems directly in . However, the indistinguishability of particles 
implies a strict symmetry requirement on which states of "Hf^ correspond to realisable physical 
states. In the case of bosons, the state must remain unchanged under an interchange of particle 
indices. For example, a two-particle state \ip) = Cij\i)\j) must have the property Cij = Cji] 
in a continuous basis this takes the more familiar form iIj{xi,X2) = 4j{x2,xi). We note that 
the correct derivation of the allowable boson and fermion symmetry rules is considerably more 
subtle than the common argument based on swapping labels in the many-particle wavefunction 
[80]. 

Taking symmetry into account, the general state for TV bosons may be written 

■ ■ • , J7v)sym = . ^^oo 5ZljPi)bP2) " " " \JPn) (2-24) 

where the sum is over all permutations P of the integers {1, . . . , A^} and rij is the number of 
particles in mode \i) as before. Such wavefunctions are usually extremely inconvenient to use 
for computational purposes due to the large number Hi^i'^iV^' unique terms. 

At this stage it is useful to introduce the number state notation. To ensure uniqueness, the 
indices in the state |ji, . . . ,jN)sym must be ordered — for example, |1, 1, 1, 2)sym and |1, 1, 2, l)sym 
are the same state due to the sum over all permutations. This redundancy can be removed by 
instead listing the occupations of all modes. We have 

|ni,ri2, . . .) = |1^_^^^^,2^_^^^^, . . .)sym, (2.25) 

ni times n2 times 

SO that, for example, |1, 1, 1, 2)sym = |3, 1, 0, 0, . . .). The state |ni, n2, . . .) has a well-defined 
number of bosons in each mode and is therefore called a number state. 



Fock space and the creation and annihilation operators 

Let J-"^ be the subspace of Hi spanned by all the correctly symmetrised number states. J-"^ 
contains all the possible states of the bosons, and therefore avoids the redundancy of non- 
physical states contained in Hi . The Fock space is constructed via the direct sum^ of all the 
J-"^ with varying A^: 

oo 

^ = 0^^. (2.26) 

N=0 



^The direct sum of two Hilbert spaces, T^i ©7^2, is the space of ordered pairs Ki/', 0)) e "Hi x ^2, with inner 
product defined by 

((^1,01)1(^2,02)) = (V'l|V'2> + (0l|02>. 
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Note the presence of the Hilbert space J-"" of zero particles which has a single basis vector 
known as the vacuum state and conventionally written |0). 

The efficiency of second quantisation arises because we can perform calculations entirely 
within J-" and the symmetry of the number states is built in from the start. The formalism 
of creation and annihilation operators is the tool used for this task. We define the creation 
operator for the ith mode by its action on the number basis: 

al\ni,...,n„...) = {rii + ly/^lm, . . . ,n, + 1, . . .) . (2.27) 

That is, the action of the creation operator aj is to create a particle in the ith mode. Note 
that this is a complete specification of dj, because we have defined its action for every number 
state and the number states are a basis for J-". From the definition we have (m + l|d^|n) = 
(n + iy^'^Smn] taking the conjugate and rearranging shows that the adjoint dj is the annihilation 
operator for mode i because it reduces the occupation by one: 

dilni, . . . ,ni, . . .) = ny^\ni, . . . , - 1, . . .). (2.28) 



A basic and important property of the creation and annihilation operators is their commu- 
tation relations. These can be computed directly from Eqs. (2.27) and (2.28) yielding 

[ai,a^]=6ij, [dj,dj]=0, [dj,dj]=0. (2.29) 

The commutation relations are vital for calculations, but may also be used as an alternate 
starting point to Eq. (2.27), resulting in the same operators. 

The operators a] are defined with respect to a particular single-particle basis {\i)}', a different 
set of operators bj. arises if we consider a different basis (Note that we are abusing the 

notation somewhat and treating kets with Greek indices as a distinct basis from those with 
Latin indices.) The change of basis law may be derived by considering the action of a] on the 
vacuum state, and using the resolution of identity 1 = in the one-particle Hilbert 

space: 

alio) = |z) = |z) = \k) = bi\0). (2.30) 

\ K / K K 

The change of basis is therefore given by 

a\ = bl and dj = 6^- (2-31) 
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Using the analogous transformation law for the continuous position basis {|x)} gives us the 
powerful concept of the boson field operator 

7^(x) = ^(x|z) ai = Y^ (2-32) 

i i 

where ipi{'x.) = (x|i) is the shape of mode i in position space. Note that we use a Greek letter 
for the field operator by convention; we could equally well have used d^(x) to emphasise the 
similarity with the operators a| for the discrete basis. We interpret the operation of the object 
'^'''(x) on a quantum state as the creation of a particle at position x. Alternatively, it may be 
thought of as a field that associates an operator on J-" to each point of space. The commutation 
relations for the continuous case are analogous to the discrete case: 

[^(x),7At(y)] = 5(x-y), [^t(x),^t(y)] = 0, [V'(x),^(y)] = 0. (2.33) 

Various expressions are shorter and more familiar when written in the position basis using ip, 
as will become clear in the following sections. 



Representation of operators on Fock space 

The creation and annihilation operators are useful because of two important properties: First, 
they act in a particularly straightforward way on the number state basis of J-", as per the defini- 
tion. Second, the operators of interest on "H^ can be written simply in terms of a] and Oj. Once 
the necessary operators are expressed in the natural basis for J-" there is no longer any need to 
consider explicitly symmetrised states defined on T-Li . This section examines the representation 
of one and two-particle operators on Fock space, largely following the development in Ref. [2]. 

We examine the case of single-particle operators first. Consider an operator q acting on "Hi. 
q may be extended to act on the space of particles as the sum 

N 

Q = J2Qi (2-34) 

i=l 

where Qi represents q operating on the ith particle^. Suppose now that q is diagonal in some 
basis {|k)} of "Hi so that Q = XIk Using the explicitly symmetrised wavefunction in 

may be written explicitly as Qi = ® g (g) i®' where i®' = i ■ ■ ■ i is the identity on HI, 

i times 

constructed from i copies of the identity operator 1 on "Hi. 
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Eq. (2.24) the action of Q on a number state may be computed. After some algebra, we find 

Q\nKi,n^2, ■ ■ ■) = (y^^qKKnn)\n,,i,n^2, ■ ■ ■) (2.35) 

K 

= l^^hi^i^ I'^'^i' "^^2, ■ ■ ■) (2.36) 

K 

^ Q = J2 = E^'^l^l'^)^-^- (2.37) 

K K 

In the second hne we have used the fact that the occupation number operator for the Kth mode 
is = Changing to the non-diagonal basis using Eq. (2.31) gives the general expression 
for any one-particle operator 

Q = Y,(i\q\jyA^a,. (2.38) 

A related procedure can be carried out for multi-particle operators, with similar results. If 
q' is a two-particle operator acting on the space "H^, an extension to l-L^ is 

N 

Q' = E Q'^^ (2.39) 

where Q'^j is q' acting on the subspace of the zth and jth particles. The operator Q' can be 
expressed in terms of the creation and annihilation operators. The procedure is similar to the 
derivation for the single-particle operators above, with some care required in handling the sum 
so that i ^ j. The general result is 

Q'=J2 ((^1 ® ® 10) ala^,akai. (2.40) 

i,j,k,l 

Note that including both orderings and (j, z) in Eq. (2.39) is redundant because g' is 

symmetric for identical particles so Q[j = Q'^^. 

The second-quantised operators we will be interested in are conveniently compact and fa- 
miliar in the position basis, as exemplified by the Hamiltonian. The single-particle part of the 
Hamiltonian^ , 



2 



^sp = -|^V^ + \/(x), (2.41) 



It is conventional to omit the hat from iJgp and other operators acting on the single-particle configuration 
space. This is somewhat inconsistent, but serves as a useful reminder of the subspace on which these operators 
act when dealing with the full Fock space. 
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extended to Fock space is 



Hi = j rfx^^(x) 



'*(x). (2,42) 



It is notable that this equation has only a single integral over the variable x rather than the 
double integral — or sum in the discrete case — which would be expected in general (cf. 
Eq. (2.38)). This is because ifgp acts only locally in the position basis, although it is not quite 
diagonal. 

If the interaction potential between a pair of particles located at x and x' is given by 
Vi^x. (x — x') , the interaction term in the TV particle Hamiltonian is 

N ^ N 

Vint (x. - X,) = 2 5Z ^i-t - ^i) • ^2.43) 

ij = l i,j=l 

This is diagonal in position space, so writing it in second-quantised form is simple: 

H2 = \jj d-K' V'nx)^nx )Mnt(x - x')^(x')^(x). (2.44) 

To summarise, the second-quantised Hamiltonian for the dilute Bose gas is given by 

H = Jd^ 7/^t(x)i/,pV;(x) + ^JJd^d^' ^l^\^)^\^')V,,,t{^ - x')^(x')^(x). (2.45) 

This Hamiltonian is the starting point for every theoretical analysis of the single species Bose 
gas [107]; it contains a complete description of the physics in the low-energy regime. 

2.2.2 The low-energy Hamiltonian 

Experiments with ultracold gases often take place in the dilute regime — the typical interatomic 
distance is much larger than the range of the potential V^nt- As a consequence, the interactions 
may be treated as two-body scattering events to a high degree of accuracy. Solving the two- 
body scattering problem involves expanding the wavefunction in a series of partial waves. 
At sufficiently low energies only the s partial wave is important and the scattering may be 
characterised by a single parameter ag called the s-wave scattering length. The scattering 
length is on the order of lOOoo for alkali atoms, where Oq is the Bohr radius, though there 
is considerable experimental flexibility in tuning the value using Feshbach resonances. For a 
detailed account of the scattering theory and related issues, we direct the reader to Ref. [96, 
Ch. 5]. 
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In practice, the dilute and low energy nature of the system implies that the theory can be 
greatly simplified by formally replacing the full interaction potential by a delta function 

T4nt(x-x') ^f/o5(x-x') (2.46) 

in the Hamiltonian, where Uq = ATrh^as/m. Strictly speaking, this replacement arises by 
integrating out the high energy states — resulting in a two-body T-matrix description — 
followed by taking the zero energy limit (see, for example, [14, §2.1]). The resulting theory 
necessarily contains a high energy cutoff that prevents the delta potential from unphysically 
scattering waves of arbitrarily high momentum. Nevertheless, we ignore this complication in 
what follows because the main approximations used later in the thesis are not sensitive to it. 
Performing the replacement of Eq. (2.46) leads to the approximate low-energy Hamiltonian: 

H = jd^ ij\^)H,J{^) + Y y rfx ^/>t(x)^t(x)V>(x)Vi(x). (2.47) 

A theory built from this Hamiltonian retains all the smooth long wavelength behaviour, while 
all short range correlations between particles are discarded. 



2.2.3 Operator equations of motion 

In section 2.2.1 we wrote operators and states for particle quantum mechanics in the con- 
venient notation of second quantisation but we said nothing about time evolution. Clearly we 
may use the Schrddinger picture in which the quantum state on Fock space evolves according 
to the usual Schrodinger equation. Equivalently we may use the Heisenberg picture of time 
evolution where the operators evolve and the quantum state is fixed. The Heisenberg equation 
of motion for an operator Ah is 

th^^^ = [AH{t),H], (2.48) 

where we have assumed that the corresponding Schrodinger picture operator As = Ah{0) is time 
independent. Note that this equation is linear — in the usual sense that linear combinations 
of solutions are new solutions — because the so-called superoperator [-,11] is linear. 

We can in principle write any desired Heisenberg picture observable in terms of the Heisen- 
berg picture field operator tpni'^jt). Therefore the full dynamics of the system are encoded in 
the behaviour of tpni'^^t), the general evolution of which is given by 

,^^h^ = ^^^(^^^t),H]. (2.49) 
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We now drop the subscripts for brevity, assuming that from now on we will work in the Heisen- 
berg picture unless otherwise specified. Using the Hamiltonian in Eq. (2.47) along with the 
commutation relations for ip"^ and ip, we obtain the operator equation of motion for ip (see, for 
example, [107]): 

,fl^}t^lll = H,J{^,t) + f/o^/'Hx,t)^/>(x,t)^(x,t). (2.50) 

The apparently nonlinear form of this equation presents a puzzle when compared to the mani- 
festly linear form of the equation of motion (2.48) for a general operator. What we have here is 
a tradeoff: We have removed the explicit dependence on the Hamiltonian operator so that the 
equation is expressed purely in terms of ip; the price is that to actually evaluate the time evo- 
lution requires computing a nonlinear function of field operators. Nevertheless, the underlying 
linearity of the solutions is preserved^ . 

The Heisenberg equations of motion arising from very simple Hamiltonians can sometimes be 
solved directly. However, including interactions as in Eq. (2.50) generally makes the problem 
analytically intractable and attacking it directly using numerics is also out of the question, 
due to the infinite number of degrees of freedom. Nevertheless, the Heisenberg equations of 
motion for the field operator give us a useful starting point for certain types of approximations, 
including the projected Gross-Pitaevskii equation discussed in section 2.4. 



2.3 The Gross-Pitaevskii equation 

The Gross-Pitaevskii equation^" (GPE) is a remarkably successful simplification of the full 
quantum field equations at the level of mean-field theory. We have already seen that for the 
ideal Bose gas all particles condense into a single quantum mode at zero temperature. The 
central approximation of the Gross-Pitaevskii theory is to assume this wavefunction is also 
reasonable for dynamics with small but nonzero interactions and excitation energies. That is, 

^For a very simple example, consider a single-mode system with Hamiltonian H = ^a^a^da representing a 
two-particle interaction. The Heisenberg equation of motion for a is 

.,da 

in— = a'aa, 
at 

which appears to be nonlinear. However, H may also be expressed in the number basis as J? = ^ ^„ n{n — 
l)|n)(n|. Writing a[t) = ^ anm{i)\n) {m\, the evolution equation for the components of a{t) is clearly linear: 

= ^ [m(TO - 1) - n{n - 1)] a„m. 

^"The GPE is commonly known as the (cubic) nonlinear Schrodingcr equation (NLSE) in other areas of 
physics and in mathematics. 
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we assume the Schrodinger picture wavefunction is 

N 

$(xi,...,x^,t) = J](^(x„t), (2.51) 

i=l 

and compute an equation of motion for (p so that the dynamics of $ approximates the full 
quantum dynamics of \1/ as closely as possible. The assumed form for the wavefunction given 
above is the time dependent Hartree-Fock ansatz. (We note that bosons in a single mode is a 
special case where the fully symmetrised Hartree-Fock ansatz is identical to the simple product 
state, also known as the Hartree ansatz.) 

There are numerous methods for deriving the equations of motion for ip. In the cold atoms 
literature the most commonly used are geared toward computing not only an equation for (/? but 
also the next order corrections, giving the so-called Hartree-Fock-Bogoliubov theories. These 
techniques may be broadly grouped into the symmetry breaking (see, for example, [24]) and 
number conserving formalisms (see, for example, [48, 88, 49]). For a discussion of the difference 
between these techniques, see Refs. [88, 107]. 

To derive only the lowest order theory, we choose a more elementary method based on a 
time dependent variational principle. 



2.3.1 The Dirac-Frenkel time dependent variational principle 

The Dirac-Frenkel time dependent variational principle [81, §11.1] has been commonly used in 
quantum chemistry and nuclear theory in the context of time dependent Hartree-Fock theory 
[70]. It is also an expedient method for deriving the time dependent GPE [96, § 7.1] and with 
that in mind we briefly outline some of its generic properties before making use of it in the next 
section. The Dirac-Frenkel variational principle is characterised by stationarity of the quantum 
effective action^^ 

r dt {^ihdt - H\^) = f dt f d^^*{ihdt- H)^ (2.52) 

J to J to J 

with respect to small variations in For an unconstrained state \& in the full Hilbert space 
T-L, the variational principle is equivalent to 

0, (2.53) 



(5* 



^"'^In path integral quantisation, one builds the quantum theory out of a classical action integrated over all 
paths. This is distinct from the full quantum "effective action" for the same system, which is a functional of 
the quantum rather than classical state. 
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where SS/S'^* is the functional derivative of the nonholomorphic function^^ 5* with respect 
to the complex variable \E'*. For a detailed discussion including how to compute with these 
derivatives in practice, see appendix A. 2. Evaluating the functional derivative gives 

§^=i'^9,-H)^, (2.54) 

which immediately yields the time dependent Schrodinger equation when set to zero. 

In contrast to the trivial case above, the variational principle is most useful when constrain- 
ing \E' to some submanifold Ai of "H; in our case M will be the set of all states of the form 
given in Eq. (2.51). In fact, the Dirac-Frenkel variational principle arises as a result of trying 
to find the "best possible" approximation of this kind: Consider some approximate evolution 
\E'(t) E Ai. According to the Schrodinger equation, the full dynamics attempts to evolve \E'(t) 
in the direction i7\E'(t), but this evolution generally takes us out of the manifold Ai and into 
the bulk of the higher dimensional space Ti. To avoid this we want an evolution direction v 
constrained to the tangent space T,i,(()A^ of Ai at \&(t) but chosen so that the error ||f — i7\I'(t)|| 
is as small as possible. Evolving \E' according to v means ihdt'^{t) = f by definition, so we 
want to minimise \\ihdt'^{t) — if\I/(t)||. This leads to the condition that the residual of the 
Schrodinger equation is perpendicular to the tangent space T^{t)M.: 

{u\ihdt - H\^!) = Q yueT^M. (2.55) 

The action in the Dirac-Frenkel variational principle follows from the orthogonality condition 
above as discussed in Ref. [81]. It is also shown that the system of equations arising from 
the variational principle automatically conserves the expected energy, E[^] = {'^\H\'^). Also 
conserved is any operator A that commutes with the Hamiltonian and keeps A'^ within T^A^: 

A^ e T^M y^eMri D{A) (2.56) 

where D{A) is the domain of A. 



2.3.2 Deriving the GPE 

Deriving the GPE is a straightforward application of the variational principle introduced in 
the previous section, with approximation manifold Ai equal to the set of all states of the form 

^^A complex-valued function of complex arguments is holomorphic if it is complex differentiable according to 
the usual limit-based definition of derivative in the complex plane. For nonholomorphic functions we must use 
a modified definition of the derivative; see section A.l. 
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given in Eq. (2.51). For convenience, we repeat the action 

'to 



S[^]= dt{'^\ihdt-H\'^). (2.57) 

Jtn 



To compute this for our particular assumed wavefunction $, we need expressions for 
and {^\ihdt\^). For the first part, we note that our wavefunction at some time t has the 
form $ = |A^, 0, . . .) in second-quantised notation, where the only occupied mode is (f{t). The 
Schrodinger picture field operators act on the state as 



iIj\N, 0, . . .) = (^(x, t)VN\N -1,0,...), (2.58) 
^"f|iV- 1,0,...) = ip*{^,t)y/N\N,0,...) (2.59) 

and therefore the expectation value of the Hamiltonian in Eq. (2.47) is 

{N,0,...\H\N,0,...) = N I rfx v?*(x,t)i/spV'(x,t) 

^ (2.60) 
+ ^N{N-1) J d^^*Mif*i^,tM^,tM^,t). 

Alternatively, this expression could be obtained by working directly with the sums over all 
particles in the first-quantised Hamiltonian^^. 

For the time derivative, it is simplest to work in the first-quantised position representation; 
we see that 

N N N 

ihdt^ = ihdtY\_v{^i^'t) = ^^^f^tV^(xj,t) ]^'y5(xj,t), (2.61) 



so that 



i=l 1=1 j=l 



N N 



{<^\zhdt\<i>) =ih I dxi--- I cixjvn<^*(xi,t)9t J]^(^(xfc,t) (2.62) 

i=l k=l 

= Nih J dxip*{x,t)dtip{x,t). (2.63) 



13 



The first-quantised Hamiltonian is 

iJ(xi,...,XAr) =^iJsp(x,) + — ^ (5(x,; -X^), 

1 — 1 

i<j 

which is equivalent to Eq. (2.47) when the number of particles is fixed at N. 
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Putting the parts of the action together, we have 

S[ip] = N j^' dt j d^^*{^,t)[imt - H,^]^{^,t) - ^{N - l)\^{^,t)\\ (2.64) 
To simphfy further, we introduce the condensate wavefunction^'^, 

0(x,t) = A^i/V(x,t), (2.65) 

so that |</)p is the particle density, and use the approximation (A^ — 1)/N ^ 1. The action is 
then given by 



S[<P] 



J^'dt J d^(t>*{^,t)[thdt-H,p](f){^,t)-^\(l){^,t)\\ (2.66) 

As a final step we apply the variational principle 6S/6(I)* = 0. Dropping the space and time 
indices for brevity, the functional derivative is 

^= [thdt-H,p]<P-Uo\M'<i^; (2.67) 

we again direct the reader to appendix A. 2 for details on computing functional derivatives. 
Setting this to zero and using the usual form for the single-particle Hamiltonian as in Eq. (2.41) 
yields the Gross-Pitaevskii equation, 

ihdtcj) = -—V^(l) + V(j) + Uo\(f)\^(l). (2.68) 
2m 



2.3.3 Validity of the GPE 

The assumptions underlying the GPE are physically reasonable at very low energies^^ and when 
the scattering length is much less than the interparticle spacing. This is borne out in the 
many successful applications to ultracold gas experiments, where the GPE is an important 
theoretical tool for understanding the dynamics of nearly pure condensates [96]. 

Nevertheless, we emphasise that the derivation given above is not particularly rigorous 
from a mathematical point of view: We have said nothing about the size of the time dependent 
errors incurred by making the product state ansatz of Eq. (2.51); neither have we shown that 

^"^Note that the condensate wavefunction is not a permissible single-particle wavefunction due to the normal- 
isation convention. 

^^It is often said that the GPE is valid at "zero temperature" . Strictly speaking, the GPE does not describe 
a state with a well-defined temperature because it is time dependent, but we may start from a true zero 
temperature state and apply a coherent excitation. 
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the product state is a good initial state for the evolution. Perhaps more severely, the effect of 
the replacement V^int(x) = ■^5(x) has also been glossed over, effectively neglecting the short 
range correlations [39] that arise from the real potential V^nt- Significant mathematical effort 
has been put into resolving these problems and convergence proofs are now available [39, 99] 
for the limit — ?■ oo with Na,^ = fixed, known as the Gross-Pitaevskii or mean-field limit. 

It is important to understand the nature of convergence in such proofs. In particular, the 
ansatz $ does not converge to the exact solution ^ of the N particle Schrodinger equation in 
the limit of large A^, 

||*(x,t) - ^(x,^)!! ^ asA^->cx). (2.69) 

Instead, it is convergence of the one body reduced density matrix^^ p(x, x') = {|\E'|^/;^(x)?/;(x') l^'^ 
to ip*{'x)(p(x'): 

||^p(x,x') -^*(x)^(x)|| -> asA^^oo (2.70) 

where ||-|| is the trace norm^''. In terms of the Penrose-Onsager criterion for condensation [95], 
this means that the system converges to a pure condensate in the mean-field limit. 

To get a feeling for the importance of using the density matrix rather than the state when 
proving convergence, consider the A^ particle state 



^' = S 



N~l 

J]<^(xfc,t) ^4^M,t)}, (2.71) 



k=l 



where ip± satisfies {(p±\ip) = and S ensures correct Bose symmetrisation of the wavefunction. 
Intuitively this state has "nearly all particles in state v?" and yet it is orthogonal to the product 
$ = '^(?^k,t) so that 11$ — $'|| is never small, even as A^ — )■ oo. On the other hand, we 
may show^*^ that the reduced density operator for $' is 

($'|^t(x)^(x')|$') = (AT - l)^*(x)(/.(x') + v^l(x)^^(x') (2.72) 

so that -^p(x, x') clearly converges to 0*(x)0(x') as A^ — )■ oo. 



^^For a well-defined number of atoms N in the pure quantum state 4*, this definition is equivalent to 
the definition p(x, x') = J • • • J (ix2 • • • dx^r \I'*(x, X2, . . . , xjv)^'(x', X2, . . . , x^v), so that the normalisation 
is / dx p(x, x) = N [96, §13.5]. This normalisation is conventional in ultracold atoms research so we use it here 
even though normalisation to the identity is more convenient for the discussion of convergence. 

i^The trace norm is defined by \\A\\ = Tr [(At^)!/^] . 

-'^^In second-quantised notation, |<i>') = |7V — 1, 1, 0, . . .) in the basis {(p, ip±, . . . }, so we have 

V;(x)|$') = {N- l)i/V(x)|iV - 2, 1, 0, . . .) + ^±(x)|iV -1,0,...), 

which allows p to be computed immediately. This is a small example of the great practical benefit of working in 
the second-quantised formalism over having to deal with explicitly symmetrised wavefunctions as in Eq. (2.71). 
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2.4 The Projected Gross-Pitaevskii equation 

The full quantum field evolution given in Eq. (2.50) contains answers to any question we might 



Unfortunately solving the operator equations of motion for i/j directly is intractable in all 
but the simplest of cases. As a consequence there is an extensive literature on methods for 
approximating the equations and/or using an alternative but more tractable form. Exact 
reformulations typically recast the problem of computing expectation values as a problem of 
stochastic sampling from some complicated probability distribution; this is the approach taken 
by the various flavours of quantum Monte Carlo and the phase space methods commonly used 
in quantum optics. Approximation methods are many and varied, but usually draw on physical 
insight about the kinds of quantum states which are of interest to reduce the complexity of a 
solution. In the following section we describe the classical field approximation in the particular 
form of the projected Gross-Pitaevskii equation (POPE), originally developed in Refs. [30, 26]. 

2.4.1 Conceptual introduction 

Under certain conditions a system governed by a quantum field theory can instead be well 
approximated using a classical field. This depends on both the state of the system and the 
kinds of measurements we wish to make. The prototypical example is the electromagnetic 
field: while the underlying theory is quantum mechanical, the classical theory — in the form 
of Maxwell's equations for the classical electric and magnetic fields — describes a huge range 
of wave phenomena very successfully. 

This begs the question: which states and types of measurements can be adequately described 
using a classical field? Qualitatively, we can say that classical fields describe only the collective 
wave-like behaviour of systems that are in reality made up of quanta, and that we should have 
an appropriately large number of those quanta so that the system is approximately continuous 
rather than discrete. More formally, this is a difficult question to answer because taking the 
classical limit of a quantum theory is a subtle business. The standard quantum states that lead 
to classical wave-like behaviour are the Glauber coherent states [52]^^. Unfortunately, there 

^^A single- mode coherent state (conventionally written \a)) is defined to be an eigenvector of the annihilation 
operator, such that a\a) = a\a) where a is a complex number. In a similar way, there is a full multimode 



ask about cold Bose gases: Given a state |\E') and a solution ipl'x^t) to the field equation, 
any observable 0{t) can be written as a function 0('?/'^(x, t), '?/'(x, t)) of the field operator; the 
expectation value of O at some time t is 




(2.73) 
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are conceptual and practical reasons that make it hard to accept coherent states as a good 
representation of the actual state found in a cold gas experiment. Conceptually, a coherent 
state is a superposition of states with different numbers of quanta, which has led to an ongoing 
debate about whether such states can even be generated in principle^" [7]. From a practical 
point of view, the number statistics of a real system is likely to include additional classical 
noise which would be inconsistent with the Poissonian number statistics of a coherent state 
[62]. These objections can be avoided if we treat the state as a classical ensemble of coherent 
states, as discussed in more detail at the end of the next section. 

The electromagnetic field is the most obvious example of a classical field arising out of 
an underlying quantum theory, but the same kind of approximation can be made for field 
theories describing ultracold gases of atoms: While the details of the equations of motion 
and the Hamiltonian differ, both theories ultimately describe the same kinds of objects — 
indistinguishable bosonic quanta having both wave-like and particle-like properties. For an 
excellent discussion of the nature of quanta, we direct the reader to Ref. [122]. 

One might worry that we now have applied the label "classical" to two mutually contradic- 
tory theories, both ostensibly describing massive particles: On the one hand, we have classical 
Newtonian particle mechanics, and on the other a classical field theory. We note that this is 
nothing other than a manifestation of wave particle duality and the fact that different classical 
limits are applicable to different situations [101]. 

To summarise, gases of cold bosons are not well described by the familiar limit of classical 
particle mechanics. As a result, we are led to the quantum theory of many particles in the form 
of a quantum field theory. Being difficult to solve, this theory is approximated — but not by 
the classical theory that we had to begin with. Instead the appropriate limit is a classical field 
theory emphasising the collective wave-like aspects of the atomic ensemble. 

2.4.2 Derivation of the PGPE 

Having said a few general things about the use of classical field theories to describe systems of 
massive particles, we turn to describing the PGPE that will be used in chapter 3 of this thesis. 
We stated earlier that classical fields are a good approximation when the number of particles 
is large. More precisely, it is the number of particles per mode (rij) which is of importance 
because the errors are of order 1/ (rij) as argued in Ref. [71]. As a result, we want each mode 
to be highly occupied, that is, (nj) 1 for alH. A system obeying this condition is necessarily 

coherent state \'^) for each complex- valued field ^(x, i) such that "01^) = ^(^j^)!^)- 

^°In fact, the possibility of creating coherent states has been called into doubt even in optics [87] where the 
photon number is not conserved. We direct the reader to Ref. [7] for a presentation of both sides of the debate 
about the reality of coherent states, and a possible resolution. 
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highly Bose degenerate (see the start of chapter 1). We note that the number of particles per 
mode is the same as the density of particles in phase space when the modes in question are the 
momentum eigenstates — that is, for a homogeneous system. For this reason the terms Bose 
degeneracy and phase space density are often used interchangeably. 

In a spatially finite system we have an infinite number of field modes, so only a finite number 
of those modes can be highly Bose degenerate as required for the classical field approximation. 
To deal with this issue, the PGPE formalism splits the set of modes into two subsets, now 
conventionally labelled the C and I regions [14]. The C or c- field region^^ is chosen to contain 
the highly occupied modes, and is simulated using a classical field. The I or incoherent region 
contains the remaining modes and is assumed to be thermalised. In thermal equilibrium the 
occupation depends on the mode energy, so the split between C and I regions is conveniently 
implemented using an appropriate energy cutoff ecut to be defined more precisely later. We 
emphasise that the choice of splitting depends on the total number of atoms, so any given 
splitting is specialised for a subset of the possible quantum states. 

There are two well known methods for deriving the PGPE. The first of these is a heuristic 
approach using what we would like to call "dequantisation" — the assertion that one may 
simply replace ipix.) by a classical field ipci'^) in the operator equations of motion. This is the 
approach taken, for example, in Ref. [30] and will be discussed further below due to its intuitive 
appeal. 

The second approach uses the truncated Wigner function formalism to derive the PGPE as 
originally described in Ref. [47] and reviewed in depth in Ref. [14]. The Wigner function method 
is attractive because it puts the classical field approximation on firmer mathematical ground and 
allows for systematic treatment of additional quantum behaviour as well as interactions with 
the I region. Nevertheless, we will not need these additional features and there is a significant 
cost in mathematical machinery, so we will not delve into the details here. We instead direct 
the reader to Ref. [14, §2] for a clear and detailed explanation. 

For a given partitioning of the modes into sets C and I, we define a pair of projection 
operators onto the subspaces spanned by the sets of mode functions {(pi : z G C} and {(pi : i El}, 
respectively: 



Note that the projectors act on spatial functions /(x), as distinct from the creation and anni- 
^"'^The C region has also been called the coherent region in earlier papers. 




(2.74) 



(2.75) 
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hilation operators that act on the Fock space. The field ip may then be spht according to 

^/'(x,t) = ^/'c(x,t) + ^/'i(x,t), (2.76) 

where ipc = Vc^^'} and ipi = Vi{ip}. At this stage no approximation has been made, and one 
can easily derive a coupled pair of equations of motion for ipc and ipi that taken together are 
equivalent to the full evolution. Applying the projector Vc to the field equation (2.50) gives 



V'lV^c^/'c + '2tphi^ii^c+ ^2 77) 



where it is assumed that the basis (pi is an eigenbasis of the single-particle Hamiltonian so that 
Vc commutes with H^p giving Vc{Hsp4j} = H^pVcW = -f^sp^c- 

We now proceed by the rather gross approximation of discarding the terms coupling the I 
region to the evolution of tpc to obtain 

ih^ = H^ptPc + UoVc{i'l4^c^Pc}. (2.78) 

This has the great advantage of providing a closed system which makes numerical work much 
simpler, but is hard to justify in general (see [26] for further comments). Nevertheless, numerical 
experience shows that the isolated C region as described by the PGPE evolves to thermal 
equihbrium [30], and a range of useful results have been obtained (see [14, §3.3-4] for several 
examples). We note that the step of discarding all interactions with the I region is a defining 
feature of the PGPE regardless of which derivation is used. Retaining these terms is necessary 
in many cases, in particular for realistic simulations of condensate formation. In such cases it 
is necessary to use a more powerful alternative, a role filled by the so-called stochastic PGPE 
(SPGPE) [47]. 

The final step in the derivation is to make the classical field approximation, that is, to take 
the classical limit of the field. On a purely formal level this is achieved by simply making the 
replacement ipc{'^,t) ipc{'^,t) in the equations of motion^^, where ipc is a complex- valued 



22 



Equivalently, this is a replacement of each mode operator with a complex amplitude Ci for all i G C. 
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classical field. We then obtain the PGPE, 

i^sp^c + f/oPcjV'cV'cV'c}, (2.79) 



which describes classical evolution of the highly occupied field modes. While this derivation is 
expedient and intuitive, it leaves much to be desired from a mathematical standpoint. 

The most basic physical justification for the replacement ipc — )■ ipc seems to be to regard 
it as a "dequantisation" or the opposite of canonical quantisation^'^. The recipe of canonical 
quantisation produces a quantum theory from a classical one by replacing the classical canonical 
position and momentum variables with operators that satisfy the usual commutation relations 
(see, for example, [41, §1.4]). We imagine dequantisation as the opposite process — forming a 
classical limit by replacing the quantum field with a classical one. 

One conceptual problem with this story is that the nature of the quantum and classical 
fields are completely different. As discussed in Ref. [122, Ch. 5], a classical field theory uses 
field equations to describe the state of the system but the quantum field is not a state-like 
object. Instead the quantum field is more akin to a Green's function: a solution ^/'(x, t) encodes 
all possible evolutions independently of any particular quantum state. In light of this, it is 
confusing to imagine replacing or somehow approximating the field operator by a classical field. 
Clearly some further mathematical justification is required. 

The classic mathematical procedure for taking the classical limit is to argue that the class of 
"classical-like" quantum states of interest are well approximated by coherent states [52]. From 
this viewpoint the replacement ipc — ^ i'c is implemented by assuming the state is coherent and 
taking the expectation value of both sides of Eq. (2.78). For example, on the left hand side we 
have 



^^9^c(x,t) 



dt 



dt at dt 



where '?/'c(x, t)|^E') = ?/'c(x, t)|\E') because l^') is a coherent state^^. This approach is known as 
the broken symmetry or mean-field approach since the average of the quantum field {ipc) = "ipc 
is nonzero for a coherent state. 

The use of coherent states is attractively simple and provides a straightforward way to take 
the classical limit. On the other hand, the coherent states are not obviously adequate for cold 



^^Canonical quantisation has a somewhat unfortunate name: It is not the canonical method of quantisation, 
but a method of quantisation making use of classical canonical coordinates. 

^^The standard notation V'c and tpc is somewhat unfortunate here, since it suggests that the state-like object 
■00 arises as an approximation to the quantum field V'c- This is not the case: it is more accurate to view ipc 
as a representation of the state l^*). A less confusing notation might be to write ^'(x, t) in place of "00 (x, t) so 
that i/ic|*) = *(x,i)|*). 
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atom experiments, for the reasons noted in Sec. 2.4.1. Luckily, we can avoid those particular 
objections if we generalise the state to a statistical mixture of coherent states as in Ref. [71]. 
The density matrix is then expressed as 



where P{a) is a distribution^^ over the set of classical fields a. With this representation at 
hand, computing the dynamics involves evolving each coherent state independently according 
to the classical field approximation [71]. Observables become ensemble averages at time t over 
the ensemble described by P. 

The procedure outlined above is not a particularly satisfying mathematical account of the 
PGPE theory, but it provides a flavour for the kinds of physical arguments that have been used 
in the field. With that in mind, we again direct the reader to Ref. [14] for a description of the 
Wigner function formalism that puts the method on more solid theoretical ground: It shows us 
more clearly which ensemble we should be using, which terms must be neglected in forming the 
equations of motion, and directions for extending the theory to deal with physical situations 
where the PGPE is not adequate. 



2.4.3 Ergodicity and thermal averages 

We are often interested in the properties of the system at a particular temperature — that 
is, the expectation values of observations of the thermal ensemble. To compute such thermal 
expectation values it is sufficient to sample any one of the standard statistical ensembles at 
the desired temperature^^. The microcanonical ensemble is simplest to deal with in this case 
— microcanonical averages should be taken over the hypersurface of constant energy, and the 
PGPE is energy conserving. If the classical energy is given by 



^^As written here, P is actually the density in the Glauber-Sudarshan P representation [52], and as a result 
is much more general than one would guess at first sight. Not all density matrices have a well behaved P 
representation (in general P can be negative and extremely singular [16]), but this is not a problem for the 
states of interest since — roughly speaking — the more classical a state is, the better behaved is the associated 
P distribution. We note that the Wigner function mentioned above is equal to the P distribution smoothed by 
a Gaussian convolution. 

^®Thc various statistical ensembles arc equivalent in the thermodynamic limit. It is worth keeping in mind 
that they have different fluctuation properties for the mesoscopic numbers of particles we are dealing with here. 




(2.81) 




(2.82) 
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then microcanonical averages are taken over all ipc with respect to the phase space density 
[123] 

{const. where HcM'c] = E and other macroscopic constraints are satisfied 
elsewhere. 

(2.83) 

We mention other macroscopic constraints because there may be additional constants of motion 
and we want any averages to take these into account. For example, the momentum is conserved 
in a homogeneous system in which case the ensemble of interest includes only microstates with 
the system at rest. 

The crucial step in computing thermal averages is to assume that the PGPE is ergodic [30] 
which allows us to convert phase space averages into time averages: 

(F) = j Vi,c P[^c; E]F[i,c] (2.84) 
= lim - / dt F[i)c{t)] (2.85) 

where F is some functional of the classical field representing an observable. In practice the time 
average is implemented numerically by sampling the motion at discrete time intervals {tj} and 
forming the sum 

mE^[^c(^^-)] (2.86) 
i=i 

for large M. This recipe provides a convenient and efficient way to sample states from the 
microcanonical ensemble. 

Not all quantities may be easily written as functionals of the field ipc. In particular we 
note that derivatives of entropy such as the temperature (T) and chemical potential (/ic) 
are calculated by time- averaging appropriate quantities constructed from the Hamiltonian in 
Eq. (2.82) using the Rugh approach [110]. The detailed implementation of the Rugh formalism 
for the PGPE is rather technical and we refer the reader to Refs. [29, 27] for additional details 
of this procedure. 

It is instructive to connect the discussion of ergodic averaging to the Wigner function version 
of the PGPE derivation. For a thermal state the initial Wigner function is very delocalised, 
in contrast to the near delta function required for a straightforward single-trajectory interpre- 
tation. Sampling directly from such a distribution is difficult except in the limit of very low 
or high temperatures where the Hamiltonian can be approximately diagonalised [14]. Happily, 
the PGPE with ergodic averaging avoids this problem by relying on the dynamics to sample 
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the distribution correctly. 



2.4.4 Treatment of the I region 

At temperatures near the BEC transition the number of atoms in the I region is significant, as 
observed in chapter 3 (see also [15] for three-dimensional systems). For realistic comparisons 
with experiment these need to be taken into account. 

A simple way to deal with the I region atoms is to assume a gas of uncorrelated semiclassical 
bosons in thermal equilibrium, interacting only indirectly via the average particle density [10, 
28]. Under these assumptions the I region can be described by an approximate single-particle 
Wigner function Wi. (Note that this is a sm(7/e-particle Wigner function, distinct from the 
multi-particle Wigner function referred to previously in section 2.4.2.) The appropriate single- 
particle Wigner function is positive and can be interpreted as a classical probability distribution 
over phase space: 



In this equation D is the number of dimensions, k is the wavevector and jj = /xc + ^UqUi is the 
chemical potential. The Hartree-Fock energy in this expression is given by 

EuFiK x) = ^ + \/(x) + 2Uo [nc(x) + ni(x)] , (2.88) 

where nc and rii are the densities of the C and I region atoms [10]; the temperature and 
chemical potential are calculated from the C region. When the potential V is not constant the 
unknown density ni is spatially varying and must be calculated self-consistently [28] using 



~'i?HF(k,x)>ecut 

This complication disappears in the homogeneous case relevant to the work in chapter 3. 

With the semiclassical Wigner function at hand, any observable -F(k, x) may be calculated 
using a phase space average 




(2.87) 





(2.89) 




(2.90) 



Note that the nontrivial region of integration satisfying £'HF(k, x) > ecut is chosen to avoid 
counting atoms that have already been taken into account via the C region simulation. 
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2.4.5 Comparison with the GPE 

We briefly contrast the PGPE formalism with the GPE, as the relationship between these is a 
recurring source of confusion: Given such apparently similar equations, why is it claimed that 
the PGPE describes all modes at temperatures up to the order of the transition temperature, 
while the GPE describes only the condensate mode at zero temperature? To further empha- 
sise the point of similarity, recall that there are various ways to derive the GPE with subtly 
different interpretations as to the exact state described. In Sec. (2.3) we presented the GPE 
as arising from the product state ansatz of Eq. (2.51), but it can just as easily be derived by 
assuming a time dependent coherent state^^. This second approach is confusingly similar to 
the "dequantisation" step in our derivation of the PGPE when interpreted in terms of coherent 
states. 

The real heart of the difference between these formalisms is the ensemble average that 
must be taken when using the PGPE to describe finite temperature states of the field. While 
this is implemented as an ergodic average, certain non-thermal states could also in principle be 
simulated by considering a full set of trajectories drawn from an initial non-thermal distribution. 
Averages would then be taken over this set of trajectories at some time^*. At finite temperature 
there must be some portion of non-condensed atoms, and indeed the PGPE is able to simulate 
these: The ensemble average results in a single particle density matrix with a largest eigenvalue 
significantly less than one; this is a non-pure condensate in the Penrose- Onsager sense. 

In contrast, the GPE simulates a single trajectory representing the motion of the condensate: 
The single particle density matrix becomes arbitrarily close to pure in the GP limit, as made 
precise in Sec. 2.3.3. The types of states that the GPE simulates are thus small perturbations 
of the zero temperature stationary solution^^. 

While the projection operator is required to provide control over the set of modes to be 
approximated, the arguments above show that it is not at the heart of the conceptual differences 
between the GPE and PGPE. Nevertheless, we emphasise that a projector is essential to clearly 

^^Thc sense in which the GPE approximates the true state of the system is addressed in Sec. 2.3.3. We point 
out here that it is consistent for both the coherent and product state ansatze to approximate the true many 
particle state in the relatively weak sense of Eq. (2.70). 

Furthermore, distinguishing between these apparently very different states in an experiment can be surpris- 
ingly difficult. This is nicely demonstrated in Ref. [125] that analyses an interference experiment, finding that 
the difference in interference fringe contrast is of order and very difficult to measure for even moderate 

numbers of particles N. This is reminiscent of the approximation N/{1 + N) « 1 that must be made in deriving 
the GPE via the product state ansatz. 

^^Indeed, this is the approach generally taken by the whole group of powerful methods arising from the 
truncated Wigner function formalism. 

^^Non-stationary solutions are not, strictly speaking, at zero temperature since simple definitions of thermal 
equilibrium imply a stationary state. A dynamical definition of temperature such as given in Ref. [110] would 
presumably assign a small but nonzero temperature to such states. 
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define the set of modes in tlie I region. Witfiout careful treatment of these high energy modes, 
one can only hope to achieve qualitative agreement with experimental observations. 



Chapter 3 

Vortex pairing in two-dimensional Bose gases 



In this chapter we investigate finite temperature 2D Bose gases using the PGPE, with a view to 
understanding the relation between BEC and BKT physics in finite-sized systems. We calculate 
several physical properties including the amount of vortex pairing, the condensate and superfluid 
fractions, and the functional form of the spatial correlations. We also relate our simulation to 
the experimental measurements described in Ref. [60]. 

3.1 Introduction 

As discussed in section 1.2, a two-dimensional homogeneous Bose gas does not undergo the 
BEC transition. Nevertheless, 2D Bose gases do display superfluid behaviour in the presence 
of interactions, due to a vortex pairing phase transition known as the Berezinskii-Kosterlitz- 
Thouless (BKT) transition. Evidence for the BKT transition has been found in experimental 
realisations of the 2D Bose gas in several studies [121, 60, 77, 117, 23]. We make particular 
note of the experiment described in Ref. [60], which was carried out at ENS in Paris, and which 
we will refer to when choosing parameters for our study. 

Experiments in the 2D regime present a new challenge for theory as strong fluctuations 
invalidate mean-fleld theories (see, for example, [105, 106, 104, 12, 83, 113, 50, 116]), and only 
recently have quantum Monte Carlo [68, 67] and classical fleld (c-fleld) [118, 119, 13] methods 
been developed that are directly applicable to the experimental regime. 

In the current chapter we study a uniform Bose gas of flnite spatial extent and parameters 
corresponding to current experiments. To analyse this system we use the PGPE, which is well 
suited to studying flnite temperature Bose flelds with many highly occupied modes. We examine 
two important applications: First, we provide a quantitative validation of the interference 
technique used in the ENS experiment to determine the nature of two-point correlation in the 
system. To do this we simulate the interference pattern generated by allowing two independent 
2D systems to expand and interfere. Applying the experimental fltting procedure to analyse the 
interference pattern, we can extract the inferred two-point correlations which we then compare 
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against the in situ correlations that we calculate directly. Second, we examine the correlations 
between vortices and antivortices in the system to directly quantify the emergence of vortex- 
antivortex pairing in the low temperature phase. A similar study was made by Giorgetti et al. 
using a semiclassical field technique [51]. We find results for vortex number and vortex pair 
distributions consistent with their results, and we show how a coarse-graining procedure can 
be used to reveal the unpaired vortices in the system. 



3.2 Formalism 



Here we consider a dilute 2D Bose gas described by the Hamiltonian 



H= /d^x^tJ + ^ I d^^^^^i^, (3.1) 



2m J 2m 

where x = (x, y). This is simply Eq. (2.47), but specialised to the two-dimensional homogeneous 
case. We take the two-dimensional geometry to be realised by tight confinement in the z 
direction that restricts atomic occupation to the lowest z mode. The dimensionless 2D coupling 
constant is 



vSyra 

9= , (3.2) 

with ttz the spatial extent of the z mode^ and a the s-wave scattering length. We will assume 
that Qz ^ a so that the scattering is approximately three-dimensional [97], a condition well- 
satisfied in the ENS and NIST experiments [121, 60, 77, 23]. For reference, the ENS experiment 
reported in Ref. [60] had g ~ 0.15, whereas in the NIST experiments g ~ 0.02 [23]. 

In contrast to experiments we focus here on the uniform case; no trapping potential in the 
xy-plane is considered. We perform finite-sized calculations corresponding to a square system 
of size L with periodic boundary conditions. Working in the finite-size regime simplifies the 
simulations and is more representative of current experiments. We note that the thermodynamic 
limit corresponds to taking L — )■ oo while keeping the density, n = (ip'^i^), constant. 



3.2.1 Review of BKT physics 

The BKT superfluid phase has several distinctive characteristics, which we briefly review. 
^For example, for tight harmonic confinement of frequency Uz we have = y/h/mujz- 
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First-order correlations 

Below the BKT transition the first-order correlations decay according to an inverse power law: 

^W(x,x') oc ||x-x'|r". (3.3) 

Systems displaying such algebraic decay are said to exhibit quasi-long-range order [19]. This 
is in contrast to both the high temperature (disordered phase) in which the correlations decay 
exponentially, and long-range ordered case of the 3D Bose gas in which g^^'^ — j- const, for 
llx — x'll — > oo. 



Superfiuid density 

Nelson and Kosterlitz [92] found that the exponent of the algebraic decay is related to the ratio 
of the superfiuid density and temperature. To within logarithmic corrections 

aiT) = ^ (3.4) 

where ps is the superfiuid density and AdB is the thermal de Broglie wavelength (Eq. (1.1)). 
Furthermore, Nelson and Kosterlitz showed that this ratio converges to a universal constant as 
the transition temperature, Tkt, is approached from below: lim^^^-^ «(^) = 1/4 (i.e., p^A^g = 
4). Thus, the superfiuid fraction undergoes a universal jump from ps(Tj|rp) = to ps(T£rp) = 
4/A^g as the temperature decreases through Tkt- 

Vortex binding transition 

Another important indicator of the BKT transition is the behaviour of topological excitations, 
which are quantised vortices and antivortices in the case of a Bose gas. A single vortex has 
energy that scales with the logarithm of the system size. At low temperatures this means 
that the free energy for a single vortex is infinite (in the thermodynamic limit), and vortices 
cannot exist in isolation. As originally argued in Ref. [75], the entropic contribution to the free 
energy also scales logarithmically with the system size, and will dominate the free energy at 
high temperatures allowing unbound vortices to proliferate. This argument provides a simple 
estimate for the BKT transition temperature. 

Although unbound vortices are thermodynamically unfavoured at T < Tkt, bound pairs 
of counter- rotating vortices may exist because the total energy of such a pair is finite^. This 
leads to a distinctive qualitative characterisation of the BKT transition: as the temperature 



^The vortex-antivortex pair energy depends on the pair size rather than the system size. 
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increases through Tkt pairs of vortices unbind. 

Location of the BKT transition in the dilute Bose gas 

While the relation Ps(^kt) ~ ^/A^g between the superfluid density and temperature at the 
transition is universal, the total density, ra, at the transition is not. General arguments [103, 
72, 43] suggest that the transition point for the dilute uniform 2D Bose gas is given by 

(nA^B)KT = In (^0 , (3.5) 

where ^ is a constant. Prokofev, Ruebenacker and Svistunov [105, 106] studied the homogeneous 
Bose gas using Monte Carlo simulations of an equivalent classical 0'* model on a lattice. Using 
an extrapolation to the infinite-sized system, they computed a value for the dimensionless 
constant, ^ = 380 ± 3. By inverting Eq. (3.5), we obtain the BKT critical temperature for the 
infinite system 

^ mkB In i^hymg)' ^^'^^ 
We use the superscript oo to indicate that this result holds in the thermodynamic limit. 



3.3 Method 

3.3.1 c-field and incoherent regions 

We briefly outline some specifics regarding how the PGPE formalism described in section 2.4 is 
applied to the two-dimensional homogeneous problem. In the homogeneous case, the fields ipc 
and ipi are defined as the low and high energy projections of the full quantum field operator, 
separated by the cutoff wave vector K. In our theory this cutoff is implemented in terms of the 
plane wave eigenstates {(pnipi.)} of the time-independent single-particle Hamiltonian, that is, 

^n(x) = -^e-'^-^ (3.7) 
kn = ^n, (3.8) 

with n = {uxyUy) G Z^. The fields are thus defined by 

V'c(x) = ^CnV3n(x), (3.9) 
neC 

V'i(x) = ^an(^n(x), (3.10) 

nel 
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where the Bose annihilation operators, the complex amplitudes, and the sets of 

quantum numbers defining the regions are 

C = {n: IlknII < K}, (3.11) 
I = {n: ||k„|| > K}. (3.12) 

Choice of C region 

In general, the apphcability of the PGPE approach to describing the finite temperature gas relies 
on an appropriate choice for K, so that the modes at the cutoff have an average occupation 
of order unity. In this work we choose an average of five or more atoms per mode using a 
procedure discussed in appendix B.l. This choice means that all the modes in C are appreciably 
occupied, justifying the classical field replacement c^- In contrast the I region contains 

many sparsely occupied modes that are particle-like and would be poorly described using a 
classical field approximation. Because our 2D system is critical over a wide temperature range, 
additional care is needed in choosing C. Typically strong fiuctuations occur in the infrared 
modes up to the energy scale fi^gn/m. Above this energy scale the modes are well described 
by mean- field theory (see, for example, the discussion in [74, 105]). For the results we present 
here, we have 

-7^ ^ n 3.13 

2m m 

for simulations around the transition region and at high temperature. At temperatures well 
below Tkt, the requirement of large modal occupation near the cutoff competes with this 
condition and we favour the former at the expense of violating Eq. (3.13). 

PGPE treatment of C region 

Specialising the PGPE (Eq. (2.79)) to 2D we have the equation of motion for ipc 

= -^^c + —Vc {l^cl ^c} , (3.14) 

where the projection operator 

Pc{F(x)} = $^(p„(x) [ dV<(x')F(x'), (3.15) 
nec ^ 

formalises our basis set restriction of ipc to the C region. The main approximation used to 
arrive at the PGPE is to neglect dynamical couplings to the incoherent region [26]. 

We assume that the evolution under Eq. (3.14) is ergodic [30], so that the microstates {ipc} 
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generated through time evolution form an unbiased sample of the equilibrium microstates. 
Time-averaging can then be used to obtain macroscopic equilibrium properties. We generate 
the time evolution by solving the PGPE with three adjustable parameters: (i) the cutoff wave 
vector, that defines the division between C and I, and hence the number of modes in the C 
region; (ii) the number of C region atoms, A^^c; (iii) the total energy of the C region, Eq. The 
last two quantities, defined as 




(3,17) 



are important because they represent constants of motion of the PGPE (Eq. (3.14)), and thus 
control the thermodynamic equilibrium state of the system. 

Obtaining equilibrium properties for the C region 

To characterise the equilibrium state in the C region it is necessary to determine the average 
density, temperature and chemical potential, which in turn allow us to characterise the I region 
(see section 3.3.2). These and other C region quantities can be computed by time-averaging as 
described in section 2.4.3. For example, the average C region density is given by 

^ Ms 

'^cW^M El^c(x,t,)l% (3.18) 

where {tj} is a set of Ms times (after the system has been allowed to relax to equilibrium) at 
which the field is sampled. We typically use 2000 samples from our simulation to perform such 
averages over a time of ~ 16 s. Another quantity of interest here is the first-order correlation 
function, which we calculate directly via the expression 




(3.19) 



The temperature (T) and chemical potential (/ic) are computed using the Rugh approach [110] 
which was briefly touched upon in section 2.4.3. 

A major extension to the formalism of the PGPE made in this thesis is the development 
of a method for extracting the superfiuid fraction, p^, from these calculations. For this we 
use linear response theory to relate the superfiuid fraction to the long wavelength limit of the 
second order momentum density correlations. An extensive discussion of this approach, and 
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the numerical methods used to implement it, are presented in chapter 4. 
3.3.2 Mean-field treatment of I region 

Occupation of the I region modes, iVj, accounts for about 25% of the total number of atoms at 
temperatures near the phase transition. We assume a time- independent state for the I region 
atoms defined by a Wigner function [91], allowing us to calculate quantities of interest by 
integrating over the above-cutoff momenta, k > K [10, 28]. 

Our assumed Wigner function corresponds to the self-consistent Hartree-Fock theory as 
applied in Ref. [28]. Specialising Eq. (2.87) to two dimensions, this is 

^^^'-'")=(2^ e(^H.(.c)-')A.T_r (3.20) 

where 

£;hf k = ^ + -{nc + ni , 3.21 

Im m 

is the Hartree-Fock energy, ni is the I region density, and /i = fic + 2h'^gni/m is the chemical 
potential (shifted by the mean- field interaction with the I region atoms). Note that the average 
densities are constant in the uniform system, so lVi(k, x) has no explicit x dependence, however, 
we include this variable for generality when defining the associated correlation function. 
The I region density appearing in Eq. (3.21) is given by 



ni= / d'kWi(k,x), (3.22) 

with corresponding atom number A^i = riiL^; total number is simply 

N = Nc + Ni. (3.23) 

An analytic expression for rii and simplified procedure for numerically calculating the first-order 
correlation function of the I region atoms, can be obtained by taking integrals over the 
phase space. These results are discussed in appendix B.2. 



3.3.3 Equilibrium configurations with fixed T and N 

Generating equilibrium classical fields with given values of Eq and Nc is straightforward since 
the PGPE simulates a microcanonical system (see appendix B.1.3). However, we wish to 
simulate systems with a given temperature and total number. As described in the preceding 
two sections these can only be determined after a simulation has been performed. In appendix 
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B.l we outline a procedure for estimating values of Ec and A'^c for desired values of and T 
based on a root finding scheme using a Hartree-Fock-Bogoliubov analysis for the initial guess. 

3.4 Results 

We choose simulation parameters in analogy with the Paris experiment of Hadzibabic et al. [60]. 
This experiment used an elongated atomic cloud of approximately 10^ ^^Rb atoms, with a spatial 
extent (Thomas-Fermi lengths) of 120 /im and 10 fim along the two loosely trapped x and y 
directions. The tight confinement in the z direction was provided by an optical lattice. 

Although our simulation is for a uniform system, we have chosen similar parameters where 
possible. Our primary simulations are for a system in a square box with L = 100 ^m, with 
4x10^ ®''Rb atoms. We also present results for systems with L = 50 fim and L = 200 fim at 
the same density in order to better understand finite-size effects. All simulations are for the 
case of g = 0.15 corresponding to the experimental parameters reported in Ref. [60]. 

The cutoff wave vector K varied with temperature to ensure appropriate occupation of the 
highest modes (see section 3.3.1). For the 100 /im system, the number of C region modes 
ranged between 559 at low temperatures to 11338 at the highest temperature studied. 

3.4.1 Simulation of expanded interference patterns between two sys- 
tems 

In order to make a direct comparison with the experimental results of Ref. [60], we have 
generated synthetic interference patterns and implemented the experimental analysis technique. 
Our simulated imaging geometry is identical to that found in Ref. [60], with expansion occurring 
in the 2;-direction. The interference pattern is formed in the xz-plane via integration of the 
density along the y-direction ("absorption imaging"). 

Our algorithm for obtaining the interference pattern due to our classical field is very similar 
to that presented in Ref. [61]. Our above-cutoff thermal cloud is taken into account separately. 
We consider a pair of fields ip''(^\x^y),ip^\x^y) from different times during the simulation, 
chosen such that the fields can be considered independent. The 3D wavefunction corresponding 
to each field is reconstructed by assuming a harmonic oscillator ground state in the tight- 
trapping direction. These two reconstructed fields are spatially separated by A = 3 fim, 
corresponding to the period of the optical lattice in Ref. [61]. 

Given this initial state, we neglect atomic interactions and only account for expansion in 
the tightly-trapped direction. This yields a simple analytical result for the full classical field 
iljc{x,y, z,t) at later times. The contribution of the above-cutoff atoms is included by an 
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incoherent addition of intensities. The result is integrated along the ^/-direction to simulate the 
effect of absorption imaging with a laser beam, that is, 

(3.24) 

i^c^ = i^ci^^ ^) + A, y, z, t). (3.25) 

Rather than integrate the full field along the y-direction, we use only a slice of length L' = 10 
fim in keeping with the experimental geometry of Ref. [60]. 

The interference patterns, ni^{x, z), generated this way contained fine spatial detail not seen 
in the experimental images. To make a more useful comparison to experiment it is necessary 
to account for the finite optical imaging resolution by applying a Gaussian convolution in the 
a;2;-plane with standard deviation 3 fim [59]. 

In accordance with the Paris experiment, we use a 22 ms expansion time to generate interfer- 
ence patterns for quantitative analysis (see section 3.4.3). To obtain characteristic interference 
images for display in Ref. [60], the experiments used a shorter 11 ms expansion [59]. We exhibit 
examples of interference patterns at various temperatures in Fig. 3.1, for this shorter expansion 
time. These images show a striking resemblance to the results presented in Ref. [60]. 



ni, 



dy 



1^, 



'^c\x,y,z,T) 



+ ni{x,y,z,T) 



3.4.2 Condensate and superfluid fractions 

For a 2D Bose gas in a box we expect a nonzero condensate fraction due to the finite spacing 
of low-energy modes. A central question is whether we can observe a distinction between the 
crossover due to Bose condensation and that due to BKT physics. To address this question we 
have computed both the condensate and superfluid fractions from our dynamical simulations. 

The condensate fraction in a homogeneous system is easily identified as the average fractional 
occupation of the lowest momentum mode. This is directly available from our simulations as a 
time average of the k = mode of the classical field, 

fc = {cM /N. (3.26) 

Superfluidity is the macroscopic tendency for some fraction fs of certain fluids to flow with- 
out apparent viscosity. The task of connecting this phenomenology to the microscopic theory 
is not trivial, so extracting the superfluid fraction from dynamical classical field simulations 
provides a more difficult challenge. Our approach relies on linear response theory to connect 
the superfluid fraction with the long wavelength limit of the second order momentum density 
correlations. The details of this technique are presented separately in chapter 4. 



48 



Chapter 3. Vortex pairing in two-dimensional Bose gases 




20 40 20 40 20 40 20 40 20 40 

X [|j.m] X [[im] X [|j.m] x [[im] x \iim] 



Figure 3.1: Synthetic interference patterns generated from the 50 /im grid by simulation of 
the experimental procedure of Ref. [60]. (a) At low temperatures, T ^ O.STkt, the interference 
fringes are straight, (b) Just below the transition temperature, T ^ 0.95Tkt, the fringes 
become wavy due to decreased spatial phase coherence. Phase dislocations become common 
at temperatures above the transition, (c) T ^ I.OSTkt, and (d) T ^ I.ITkt- These "zipper 
patterns" indicate the presence of free vortices, (e) When simulation of the finite imaging 
resolution is disabled, the zipper patterns from the field in sub figure (d) are no longer clearly 
visible; the high frequency details obscure the phase information without providing obvious 
additional information about the existence of vortex pairs. 

Figure 3.2 compares the results for the superfluid and condensate fractions computed on 
the 100 /im grid. These results are qualitatively similar to the results for the larger and smaller 
grids. In particular, we note that there is no apparent separation between temperatures at which 
the superfluid and condensate fractions fall to zero. Also shown in Fig. 3.2 is the condensate 
fraction for the ideal Bose gas confined to an identical finite-size box in the grand canonical 
ensemble. The large shift between ideal and computed transition temperatures indicates the 
effect of interactions in the 2D system. Because the average system density is uniform, this 
large shift is to due to critical fluctuations (also see [74]). 

In our calculations we identify the transition temperature, Tkt, as where the superfluid 
fraction falls off most rapidly (i.e., the location of steepest slope on the fg versus T graph; see 
Fig. 3.2). As the system size increases, this transition temperature moves toward the value for 
an infinite-sized system, [105]. This effect is illustrated by the behaviour of the superfluid 
fraction in Fig. 3.3. 

Also shown in Fig. 3.3 is an alternative calculation of the superfluid fraction based on 
Eq. (3.4). The two methods are expected to match for temperatures at and slightly below the 
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Figure 3.2: Condensate fraction (solid dots) and superfluid fraction (crosses) as functions of 
temperature for the 100 /im^ grid. The transition temperature in the thermodynamic limit, 
[105], is shown as a vertical dot-dashed line. The vertical dashed line shows our estimate for the 
transition temperature in the finite system. The thick solid line is the condensate fraction for 
an ideal Bose gas in the grand canonical ensemble with the same number of atoms and periodic 
spatial domain. The superfluid fraction becomes negative in places because the extrapolation 
of the momentum correlations to k = is sensitive to statistical noise at high temperature (see 
section 4.3 for details). 



transition temperature [92]. 



3.4.3 First-order correlations — algebraic decay 

Algebraic decay of the first-order correlations, as described by Eq. (3.3), is a characteristic 
feature of the BKT phase. Above the BKT transition, the first-order correlations should revert 
to the exponential decay expected in a disordered phase. 

The normalised first-order correlation function, ^(1) is defined by 



G(i)(x, X' 



'n xm xM 



where G^^^'XjX') = ^^/''''(x)^/'(x')) is the unnormahsed first-order correlation function [91]. In a 
homogeneous isotropic system g^^'' depends only the distance ||x — x'|| and we may characterise 
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Figure 3.3: Detail of the superfluid fraction near the transition temperature. Solid dots 
represent the calculation based on momentum correlations as described in chapter 4. Results 
for the largest and smallest grids are shown (left and right, respectively). The data for the 100 
/xm grid is omitted for clarity, but lies between the curves shown as expected. Open circles 
represent the calculation of the superfluid fraction from the associated fitted values for the decay 
coefficient a, via Eq. (3.4). The open circles terminate where the power law fitting procedure 
fails. 

the first-order correlations by a function of one variable, g^^\x) = g^^\\\x — x'||) = ^''^^^(x, x'). 
Direct calculation of g^^^ 

In the PGPE formalism the C and I contributions to the correlation function are additive [10], 
that is, 

G'(^)(x,x') = G'g^(x,x') + G['^(x,x'), (3.28) 

where G^^ and are defined in Eqs. (3.19) and (B.16), respectively. It is interesting to note 
that and G^^ individually display an oscillatory decay behaviour — originating from the 
cutoff — an effect which correctly cancels when the two are added together. 

Having calculated g^^\ we obtain the coefficient a by fitting the algebraic decay law, 
Eq. (3.3), using nonlinear least squares; sample fits are shown in Fig. 3.4. The fit is con- 
ducted over the region between 10 and 40 de Broglie wavelengths. The short length scale cutoff 
is to avoid the contribution of the non-universal normal atoms, for which the thermal de Broglie 
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wavelength sets the appropriate decay length. The long distance cutoff is chosen to be small 
compared to the length scale L, to avoid the effect of periodic boundary conditions on the long 
range correlations. 

The quality of the fitting procedure, and the breakdown of the expression in Eq. (3.3) at 
the BKT transition can be observed by adding an additional degree of freedom to the fitting 
function. In particular, at each temperature we fit the quadratic \io.{g^^^) = A~a\n{x) + 5\n^{x) 
and extract the parameter 5 (a ~ a is discarded). The abrupt failure of the fits can be observed 
in the inset of Fig. 3.5 as a sudden increase in the value of \S{T) \ — an effect which is in excellent 
agreement with the value of Tkt as estimated from the superfluid fraction. 




X [|a,m] 

Figure 3.4: Sample fits to the algebraic decay of g^^^ at multiples T ^ 0.77,0.93,1.01 and 
1.12 of the transition temperature. High temperatures correspond to curves at the bottom of 
the figure which have rapid falloff of g^^^ with distance. Fits are shown on a log-log scale in 
the inset to emphasise the failure of a power law in describing the behaviour of g^^^ at high 
temperature. 

Calculation of g^^'' via interference patterns 

So far a direct probe of the in situ spatial correlations has not been possible, although impor- 
tant progress has been made by the NIST group [23]. In the experiments of Hadzibabic et al. 
[60] a scheme proposed by Polkovnikov et al. [102] was used to infer these correlations from the 
"waviness" of interference patterns produced by pair of quasi-2D systems (see section 3.4.1). 
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Figure 3.5: Comparison of two methods for determining the algebraic decay coefficient a{T) 
for the first-order correlation function (^'^^•'(x, x'). The line with circle markers represents direct 
fits to g^^\ These fits fail at the transition temperature as shown by the sharply diverging value 
of |5(T)| in the inset. The filled points represent the values a'{T) obtained from a simulation of 
the experimental analysis procedure of Ref. [60], described in section 3.4.3. Horizontal dotted 
lines at 0.25 and 0.5 correspond to the expected values of a' just below and above the transition, 
respectively [60]. The vertical line is the BKT transition temperature, as estimated from the 
superfiuid fraction calculated in section 3.4.2. 



In this section we simulate the experimental data analysis method, and compare inferred pre- 
dictions for the correlation function against those we can directly calculate. This allows us to 
characterise the errors associated with this technique arising from finite-size effects and finite 
expansion time. 

To make this analysis we follow the procedure outlined in Ref. [60]. We fit our numerically 
generated interference patterns (see section 3.4.1) to the function 



F{x,z) = G{z) 



1 + c(x) COS 



27rz 

IT 



d{x] 



(3.29) 



where G{z) is a Gaussian envelope in the z-direction, c{x) is the interference fringe contrast, 
D is the fringe spacing and 9{x) is the phase of the interference pattern in the z-direction. 
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Defining the function 

C{L^) = ^I'^ dx c{x)e''^^''\ (3.30) 

tlie nature of spatial correlations is then revealed by the manner in which <||C(L^.)p) decays 
with Lx- In particular, we identify the parameter a', defined by (|C(L^.)p) oc L~'^°'' [102]. For 
an infinite 2D system in the superfluid regime (T < T^) a' = a (i.e., a' corresponds to the 
algebraic decay of correlations). For T > T^, where correlations decay exponentially, a' is 
equal to 0.5. 

Fitting <(|C(La.)p) to the algebraic decay law AL~'^"' we can determine a' . A comparison 
between a' inferred from the interference pattern and a obtained directly from is shown in 
Fig. 3.5. Both methods give broadly consistent predictions for a when T < Tkt, however our 
results show that there is a clear quantitative difference between the two schemes, and that a' 
underestimates the coefficient of algebraic decay in the system (i.e., using a' in Eq. (3.4) would 
overestimate the superfluid density). Near and above the transition temperature, where the 
fits to g^^^ fail, we observe that a' converges toward 0.5. The agreement between a and a' in 
the low temperature region improves as the size of the grid is increased. 

3.4.4 Vortices and pairing 

The simplest description of the BKT transition is that it occurs as a result of vortex pair 
unbinding: At T < Tkt vortices only exist in pairs of opposite circulation, which unbind at the 
transition point to produce free vortices that destroy the superfluidity of the system. However, 
to date there are no direct experimental observations of this scenario, and theoretical studies 
of 2D Rose gases have been limited to qualitative inspection of the vortex distributions. In 
the c-field approach vortices and their dynamics are clearly revealed, unlike other ensemble- 
based simulation techniques such as quantum Monte Carlo where the vortices are obscured by 
averaging^. This gives us a unique opportunity to investigate the role of vortices and pairing 
in a dilute Bose gas. 

We detect vortices in the c-field microstates by analysing the phase profile of the instanta- 
neous field (see appendix B.3). An example of a phase profile of a field for T < Tkt is shown 
in Fig. 3.6(a). The vortex locations reveal a pairing character, that is, the close proximity 
of pairs of positive (clockwise) and negative (counterclockwise) vortices relative to the aver- 
age vortex separation. An important qualitative feature of our observed vortex distributions 
is that at high temperatures, pairing does not disappear from the system entirely. Indeed, 

■^For example, [42] calculates the vortex density, but only indirectly via a relation with the quasiparticle 
density. 
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most vortices at high temperature could be considered paired or grouped in some manner, as 
shown in Fig. 3.6(b). Perhaps this is not surprising, since positive and negative vortices have 
a logarithmic attraction, and we observe them to create and annihilate readily in the c-field 
dynamics. However, this does indicate that the use of pairing to locate the transition may be 
ambiguous, and we examine this aspect further below. 

It is also of interest to measure the number of vortices, A^^,, present in the system as a 
function of temperature (see Fig. 3.7). At the lowest temperatures the system is in an ordered 
state, and the energetic cost of having a vortex is prohibitive. As the temperature increases 
there is a rapid growth of vortex population leading up to the transition point followed by linear 
growth above Tkt- 



Radial vortex density 

The most obvious way to characterise vortex pairing is by defining a pair distribution function 
for vortices of opposite sign. Adopting the notation of Ref. [51], this is 

Ga(r) = (p,,+ (0)p,,_(r)), (3.31) 

where p^^_|_ is the vortex density function which consists of a sum of delta spikes, 

p.,4r) = 5^ 5(r - r+) (3.32) 

for positive vortices at positions {r^}- We use the analogous definition for p^^^. The associated 
dimensionless two-vortex correlation function is 

(2) / X _ ^l.il'^) /o qq\ 

(Pt',+ (0))(Pt,,-(r)) 

(2) 

The angular average of g^J^ can be calculated directly from the detected vortex positions using 
a binning procedure on the pairwise distances ||r^ — rj||, and is shown in Fig. 3.8. 

These results quantify the effect discussed earlier: Positive and negative vortices show a 
pairing correlation that does not disappear above Tkt- The characteristic size of this correlation, 
given by twice the width of the peak feature in Fig. 3.8, is /cor ~ 3/im (taking full width half 
maximum) . 

The shape of our pairing peak is qualitatively similar to that described in Ref. [51]. However, 
in contrast to their results the width does not appear to change appreciably with temperature. 
Additional simulations show that increasing the interaction strength causes the peak to become 
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Figure 3.6: Phase profile of a c-field witli vortices indicated. Vortices witli clockwise (wliite +) 
and anticlockwise (black o) circulation. The phase of the classical field is indicated by shading 
the background between dark blue (phase 0) and light yellow (phase 27r). (a) Distinctive pairing 
below the transition at T = 207nK ^ 0.93Tkt- (b) A "vortex plasma" above the transition at 
T = 238nK ^ I.OTTkt- 
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Figure 3.7: Total number of vortices (dots) and number of unpaired vortices (circles) as a 
function of temperature near the transition. While at the transition temperature is already 
very high, becomes nonzero only close to the transition, providing clear evidence of vortex 
unbinding at work. The inset shows the variation in the total number over the full temperature 
range of the simulations. Above the transition temperature the growth in the number of vortices 
becomes linear with temperature. 

squarer and wider. It is clear that while the pair size and strength revealed in gl;±{r) does not 
change appreciably as the transition is crossed, the amount of pairing relative to the background 
uncorrelated vortices changes considerably. This background of uncorrelated vortices is given 
by the horizontal plateau gy^±{r) — )■ 1 at large r as shown in the inset. 

Revealing unpaired vortices with coarse-graining 

(2) 

The function Gl 4 (r) clearly indicates the existence of vortex pairing in the system. However, it 
does not provide a convenient way to locate the transition temperature, because a large amount 
of pairing exists both below and above the transition: The expected number of neighbours for 
any given vortex — roughly, the area of the pairing peak of (nt,^_|_)G[,^j_(r) shown in Fig. 3.8 — 
does not change dramatically across the transition. (n^,+) = (n„)/2 is the expected density of 
positive vortices. 

We desire a quantitative observation of vortex unbinding at the transition and have therefore 
investigated several measures of vortex pairing^. However, measures based directly on the full 

*For example, the Hausdorff distance (see, for example, [94, p. 105]) between the set {r^} of positive vortices 
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Figure 3.8: Angular average of the two-vortex pair distribution functions for vortices of 
opposite sign. Three temperatures centred about the transition are shown: dot markers T = 
194nK ^ 0.9Tkt, fc = 0.34; circle markers T = 217nK ^ I.OITkt, fc = 0.076; cross markers 
T = 236nK ^ I.ITkt, fc = 0.006. The vertical dotted line shows the value of the healing length 
at T = 0. The main plot shows gl^l. normalised by the positive vortex density; comparable 
magnitudes for the peaks near r = show that vortex pairing remains important over the range 
of temperatures studied, not only below the transition. The inset shows 4 in the natural 

(2) 

dimensionless units for which gl,^±{r) — )■ 1 as r — )■ oo. 

set of vortex positions seem to suffer from the proliferation of vortices at high temperature 
— an effect that tends to wash out clear signs of vortex unbinding. With this in mind, we 
have developed a procedure for measuring the number of unpaired vortices in our simulations, 
starting from the classical field rather than the full set of vortex positions. 

The basis of our approach for detecting unpairing is to coarse-grain the classical field by 
convolution with a Gaussian filter of spatial width (standard deviation) a/. This removes 
all vortex pairs on length scales smaller than o"/. Figure 3.9 shows the count of remaining 
vortices as a function of filter width, along with some examples of coarse-grained fields. For 
(7/ > /cor, the number of remaining vortices levels off and only decreases slowly with increasing 
(T/. Ultimately the number of remaining vortices goes to zero as aj — L. 

Setting the filter width to be larger than the characteristic pairing distance, /cor, yields a 



and the set {r^ } of negative vortices. 
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coarse-grained field from wliicli tlie pairs liave been removed, but unpaired vortices remain. In 
our simulations we liave /cor ~ 3 fim; we take the vortices tliat remain after coarse-graining 
with a Gaussian of standard deviation ctj = 5 /im to give an estimate of the number of unpaired 
vortices, N^. Figure 3.7 shows that A^^ becomes nonzero only near the transition, in contrast 
to A^^ which is nonzero well below Tkt- The sharp increase in at Tkt is a quantitative 
demonstration of vortex unbinding at work. 



100 




filter width [^im] 

Figure 3.9: The coarse-graining procedure: number of vortices as a function of filter width 
for a temperature near the transition. The smooth curve is an average over many realisations 
of the field, whereas the stepped curve shows typical behaviour of the number for a single 
field. Insets show the coarse-grained fields for various filter widths; the transformation removes 
vortex-antivortex pairs that are separated by approximately less than the standard deviation 
of the filter. In this example Nu = 4 unpaired vortices remain at cx/ = 5 /im. 

In the experiment of Ref. [60], the fraction of interference patterns with dislocations (see, 
for example. Figs. 3.1(c) and (d)) was measured. While isolated vortices are clearly identified 
by interference pattern dislocations, a lack of spatial resolution in experiments means that 
this type of detection method obscures the observation of tightly bound vortex pairs. The 
experimental resolution of 3 fim is broadly consistent with the scale of the coarse-graining filter 
(i.e., cTf = 5 ^m). With this in mind, we introduce the quantity Pu(T), defined as the probability 
of observing an unpaired vortex in a 50 x 50 fim control volume at a given temperature^. For 
^We choose a fixed control volume with L = 50 /im in order to compare results between simulations with 
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the 50 /iin grid we have simply Pu{T) = Pt{Nu > 1). 

Computing Pu{T) from our simulations yields the results shown in Fig. 3.10. Our results 
show a dramatic jump in pu at a temperature that is consistent with the transition temperature 
Tkt determined from the superfluid fraction calculation presented in section 3.4.2. 



3 

Q. 
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Figure 3.10: Comparison of vortex unpairing measures. The dots are our pairing measure 
based on coarse-graining the field. Circles represent the pairing as determined by the number 
of dislocations in the simulated interference patterns. This was the same method used in the 
experimental analysis of Ref. [60] and coincides remarkably well with our coarse-graining based 
measure. Both curves are consistent with the vertical line showing the transition temperature 
Tkt as determined from the superfluid fraction calculation in section 3.4.2. The inset shows 
the calculated coarse-grained pairing measure for all three grid sizes, along with vertical lines 
showing the estimates for Tkt derived from the superfluid fraction calculations. 

From the definition, we expect that p^ should be close to the experimentally measured 
frequency of dislocations. To demonstrate this relationship, we have simulated interference 
patterns (as described in section 3.4.1) and detected dislocations using the experimental pro- 
cedure of Ref. [60]: A phase gradient dO/dx was considered to mark a dislocation whenever 
\d9/dx\ > IT /A rad//im. From this we can compute the probability of detecting at least one 
dislocation as a function of temperature. As shown in Fig. 3.10, the results of this procedure 
compare very favourably with our measure of pairing based on pu- We note that inhomoge- 

diffcrcnt grid sizes. 
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neous effects in experiments probably broaden the jump in pu appreciably compared to our 
homogeneous results. 

3.5 Conclusion 

In this chapter we have used c-field simulations of a finite-sized homogeneous system in order to 
investigate the physics of the 2D Bose gas in a regime corresponding to current experiments. We 
have directly computed the condensate and superfiuid fractions as a function of temperature, 
and made comparisons to the superfiuid fraction inferred both from the first-order correlation 
function, and the experimental interference scheme. Our results for these quantities provide a 
quantitative test of the interference scheme for a finite system. 

An intriguing possibility is the direct experimental observation of vortex-antivortex pairs, 
their distribution in the system, and hence a quantitative measurement of their unbinding at the 
BKT transition. We have calculated the vortex correlation function across the transition and 
provided a coarse-graining scheme for distinguishing unpaired vortices. These results suggest 
that the dislocations observed in experiments, due to limited optical resolution, provide an 
accurate measure of the unpaired vortex population and accordingly are a strong indicator of 
the BKT transition. 

We briefly discuss the effect that harmonic confinement (present in experiments) would 
have on our predictions. The spatial inhomogeneity will cause the superfiuid transition to be 
gradual, occurring first at the trap centre where the density is highest, in contrast to our results 
where the transition occurs in the bulk. 

Bisset et al. [13] used an extension of the c-field method for the trapped 2D gas to examine 
g^^'^ and found similar results for the onset of algebraic decay of correlations at the transition. 
Their analysis was restricted to the small region near the trap centre where the density is 
approximately constant; we expect the results of our vortex correlation function and the coarse- 
graining scheme should similarly be applicable to the trapped system in the central region. 
Except in very weak traps, the size of this region is relatively small and will likely prove 
challenging to measure experimentally. 

Our results for the homogeneous gas emphasise the clarity with which ab initio theoreti- 
cal methods can calculate quantities directly observable in experiments, such as interference 
patterns. This should allow direct comparisons with experiments, providing stringent tests of 
many-body theory. 



Chapter 4 

Superfluid fraction and the PGPE 



In this chapter we describe a method for calculating the superfluid fraction from a PGPE clas- 
sical field simulation. We first present the analytical derivation, showing how the superfluid 
density arises in certain limits of the momentum density autocorrelations. This is followed by 
discussion of a numerical implementation relevant to the 2D simulations of chapter 3. 

4.1 Introduction 

Superfluidity is a famous example of macroscopic quantum behaviour, and is typically discussed 
in macroscopic terms. In particular, one characterises a superfluid by its zero viscosity; the 
ability to "flow without friction" through a narrow channel. Although such macroscopic ideas 
are easily expressed, it is not trivial to connect them to microscopic theories such as the PGPE 
formalism in an efficient way. In the following we briefly provide some relevant background 
before following with the details of our derivation in the next section. For further background 
theory we refer the reader to chapter 6 of Ref. [100] which provides an accessible overview of 
superfluid theory as relevant to experiments on ultracold Bosons. 

According to Landau's phenomenological model of superfluidity [78], it is possible to model 
a superfluid system as a "mixture" of two liquids: a superfluid part without viscosity, and a 
normal part. This idea was introduced to describe the residual viscosity which remains in liquid 
helium, even below the superfluid transition temperature. One deflnes a superfluid fraction as 
the ratio fs = Ps/p of the superfluid density ps to the total density p of the system^. 

The two-fluid model was put on flrmer ground by Putterman and Roberts in Ref. [108]. 
Starting from the equations for a single nonlinear classical fluid, they considered the presence 
of small amplitude excitations on top of a background fluid. Using only a separation of scales 
argument, they have derived kinetic equations for these thermal excitations. In the hydrody- 
namic (collision dominated) regime, the model then reduces to the Landau two-fluid model 

^Note that we follow convention and use the mass density p in the eurrent chapter rather than the number 
density n that is used in the rest of the thesis. 
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of superfluidity. The connection between this model and the classical field methods is further 
described in [112]. 

In Landau's theory, superfluidity may be predicted from the form of the energy spectrum 
of elementary excitations. Let e(p) be the energy of an excitation with momentum p measured 
with respect to a stationary background of fluid. To understand the origin of superfluidity, 
we imagine that the fluid occupies a narrow channel with walls moving at velocity u with 
respect to the fluid. In the frame of reference where the walls are stationary, the energy of 
the excitation is found to be e'(p) = e(p) — u ■ p after applying a Galilean transformation to 
the Hamiltonian^. This is negative for sufficiently large u, making the formation of excitations 
energetically favourable. However, if u is small enough, it may be the case that e'(p) > for all 
p and excitations are energetically forbidden. This leads to Landau's criterion for superfluidity: 
||u|| must be smaller than the critical velocity. 




where p = ||p||. 

It is worth noting that any system with quadratic dispersion relation e(p) oc p"^ for small p 
has a critical velocity of zero and cannot be a superfluid. In particular, this includes the ideal 
gas, where the elementary excitations are simply particles and we have e(p) = p^ /2m. On the 
other hand, introducing interactions as in the Hamiltonian Eq. (2.47) modifies the energy so 
that e(p) oc p at small p, which allows the system to support superfiuidity. For sufficiently 
weak interactions, Hamiltonian (2.47) may be approximately diagonalised via the Bogoliubov 
transformation (see, for example, [100, §4.3]), yielding the classic dispersion relation 



for the energies of elementary excitations, known as Bogoliubov quasiparticles. 

The approximate diagonalisation discards terms corresponding to the interaction of quasi- 
particles, which is a good approximation for sufficiently weak interactions and low temperatures. 
Further, this leads to a well known method for computing the superfluid fraction. At nonzero 
temperature, the Bose distribution gives the number of non-interacting quasiparticles at each 
energy. 




Vc = mm 
p 



(4.1) 



p 




(4.2) 




(4.3) 



^Note that in our notation u is the velocity of the walls with respect to the fluid, so the superfluid velocity 
is — u with respect to the walls. This is opposite from the convention used in Ref. [100]. 
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This may be used to evaluate the expected momentum density (p)^, which we attribute to the 
motion of the normal fraction, so that 



where we have used the phase space volume appropriate to three dimensions. For consis- 
tency with the next section, we use the notation to mean an expectation value with respect 
to the system where the walls are moving. After some manipulation, we arrive at an expression 
for the density of the normal fraction [100] 



This method is valid when the physical picture of non-interacting quasiparticles is valid 
— in particular, one requires weak interactions and low temperatures. In the context of the 
classical field method for ultracold Bose gases, the method has previously been used to compute 
the superfiuid fraction, see for example Ref. [127]. Unfortunately it is not valid for the system 
considered in chapter 3 for two reasons. First, we wish to compute the superfiuid fraction over 
a wide range of temperatures, from zero to slightly above the transition and the underlying 
assumptions are invalid near the transition [100, p. 66]. Second, the weakly interacting limit is 
especially difficult to reach in two dimensions because it requires the inequality lnln(l/na^) <^ 1 
to be satisfied as discussed in Ref. [43] (see also Ref. [106]). 

Given our two-dimensional system at moderate temperatures, we require a non-perturbative 
alternative to Eq (4.5). Starting from the macroscopic definition of superfluidity, one naturally 
imagines performing a time-varying numerical experiment in order to determine the superfluid 
density within the PGPE formalism. For example, in the homogeneous case we might examine 
the drag force produced when a perturbing potential is moved across the system; zero drag 
would imply a superfluid fraction of 100%. Dynamical PGPE simulations were used in a 
similar way in Ref. [119] to provide evidence for superfluidity in a 2D trapped Bose gas, by 
analysing the dynamics of the "scissors mode" oscillation. 

Deducing superfluidity from dynamical simulations is certainly possible, but is far from 
ideal: For the work presented in chapter 3 it would mean performing an additional set of 
numerically expensive simulations at every temperature of interest. Our method address this 
problem by using linear response theory to relate the superfluid fraction to the long wavelength 
limit of the second order momentum density correlations. The method is attractive because 
the momentum correlations may be extracted directly from PGPE simulations at thermal equi- 
librium. This allows the superfluid fraction to be computed from the same set of simulations 




(4.4) 




(4.5) 
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as the temperature, chemical potential, and other thermal properties; there is no need to per- 
form an expensive special purpose simulation for the sole purpose of calculating the superfluid 
fraction. 



4.2 Momentum density correlations and the superfluid 
fraction 

Our derivation is based on the microscopic theory presented in Ref. [45, p. 214], (see also [8] 
and [100, p. 96]). The central idea is to establish a relationship between (i) the autocorrelations 
of the momentum density in the simulated ensemble and (ii) the linear response of the fluid to 
slowly moving solid boundaries; (i) is a quantity we can calculate, while (ii) is related to the 
basic properties of a superfluid via a simple thought experiment. 

To connect the macroscopic, phenomenological description of superfluidity with our mi- 
croscopic theory, we make use of the standard thought experiment shown schematically in 
Fig. 4.1(b): Consider an infinitely long box, B containing superfluid, and accelerate the box 
along its long axis until it reaches a small velocity u. Due to viscous interactions with the 
walls, such a box filled with a normal fluid should have a momentum density at equilibrium 
of (p)u = pu. As above, the notation (•)u denotes an expectation value in the ensemble with 
walls moving with velocity u. 

Because the superfluid part is nonviscous, the observed value for the momentum density 
in a superfluid is less than the value pu expected for a classical fluid. As above, we attribute 
the observed momentum density, p^u, to the "normal fraction" where p„ is the normal fluid 
density. The superfluid fraction remains stationary in the lab frame, even at equilibrium and 
makes up the remaining mass with density ps = p — pn- 

In order to apply the usual procedures of statistical mechanics to the thought experiment, we 
consider two frames: the "lab frame" in which the walls move with velocity u in the x-direction 
and the "wall frame" in which the walls are at rest. 

Assuming that the fluid is in thermal equilibrium with the walls, the density matrix in the 
grand canonical ensemble is given by the usual expression p = e~^^^^~^^^ / Tr (^e~^^^'^~^^'^^ 
where is the Hamiltonian of the system in the wall frame and f3 = l/ksT. A Galilean 
transformation relates to the Hamiltonian in the lab frame, = H — u ■ P + ^Mu^, 
where P = d^x p(x) is the total momentum, M = mN is the total mass and p(x) is the 
momentum density operator at point x. The expectation value for the momentum density in 
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Figure 4.1: Thought experiment used in deriving the superfluid density. The walls move with 
velocity u in the x direction. To begin with, we imagine that the superfluid sits in a box of 
dimensions shown in (a). We later take the limit as the box walls recede to inflnity 

to get the thermodynamic limit (d). The order of the limits is critically important: the path 
(b) leads to superflow while the path (c) results in the entire fluid moving along with the walls. 



the presence of moving walls is then given by the expression 

(p(x))^ = Tr[pp(x)], 

Expanding this expression to flrst order in u yields 

(P(x)). = (P(x)) + /3((p(x)u . P) - (p(x)) (u ■ P)) , 



(4.6) 
(4.7) 



(4.^ 



where all the expectation values on the right hand side are now taken in the equilibrium ensemble 
with the walls at rest. Since (p(x)) = in our equilibrium ensemble, this simplifies to 



(p(x)), = /3(p(x)P> ■ u, 

= /3 / rfV (p(x)p(x')) ■ u, 
Jb 



(4.9) 
(4.10) 
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where p(x)p(x') is a rank-two tensor; the outer product of p(x) and p(x'). [In two dimensions 
this means (p(x)p(x')) is a 2x2 matrix for each pair of coordinates (x, x').] 



To make further progress, we wish to take the hmit as the system gets very large (we 
will notate this limit as B ^ oo). To this end, we first consider some properties of the 
correlation functions in the infinite system. The infinite system is homogeneous, which implies 
that (p(x)p(x'))^ = (p(x + r)p(x' + r))^ for any r, where (■)^ indicates an average in the 
infinite system. As a consequence, we may express the correlations — in the infinite system — 
in terms of the Fourier transform in the relative coordinate x' — x: 

(p(x)p(xO)^ = (p(O)p(x' - x))^ (4.11) 
= (2^/ ^'ke^^-(-'-'')x(k), (4.12) 

where all the important features of the correlations are now captured by the tensor 

X(k) = I dhe~^^-{mPir))^. (4.13) 

Because of the isotropy of the fluid in the infinite system, x(k) obeys the transformation law 
x(Ok) = 0~^x(k)0, for any 2x2 orthogonal matrix O. This implies that x may be decomposed 
into the sum of longitudinal and transverse parts: 

X(k) = kkxiik) + (J - kk)xt{k) (4.14) 

where k = k/k, k = ||k||, I is the identity and the juxtaposition of vectors kk represents the 
outer product as above. The transverse and longitudinal functions Xt and xi ^i-re scalars that 
depend only on the length k. 



We now return our attention to the finite system. If the finite box B is large then the 
momentum correlations in the bulk will be very similar to the values for the infinite system. 
Therefore, when x and x' are far from the boundaries, we may approximate 

(p(x)p(x')) ^ (P(x)p(x'))^ (4.15) 
= (2^/ d'^^'^^'-'-'-hi^) (4.16) 
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which in combination with Eq. (4.10) yields 

(p(x)).^/3y^rfV^ I d'ke^^<-'-h{k).n (4.17) 
= (3 J rf'k AB(k)e-*'^-"x(k) ■ u. (4.18) 

Here we have defined the nascent delta function A5(k) = d'^x' e^^'^' which has the 

property AsiX) — ?■ (5(k) as 5 — )■ oo. 

We are now in a position to carry out the limiting procedure to increase the box size to 
infinity. However, care must be taken because the simple expression limB_5.oo (p(x))^ is not well 
defined without further qualification of the limiting process i? — t- oo. 

To resolve this subtlety we must insert a final vital piece of physical reasoning. Let us 
assume for simplicity that u is directed along the x-direction, and the box B is aligned with 
the X and y axes with dimensions L^xLy. As shown in Fig. 4.1, there are two possibilities for 
taking the limits, representing different physical situations. 

On the one hand [Fig. 4.1(b)], we may take the limit oo first, which gives us an 

infinitely long channel in which superfluid can remain stationary while only the normal fraction 
moves with the walls in the x-direction. We have 

Pnii = lim lim (p(x))^ (4.19) 

Ly—^oo Lx—^oo 

= lim lim (3 /d2kAB(k)e-^''-''x(k) -u (4.20) 
= (3 lim lim x(k) ■ u (4.21) 

where we use the fact that Asi^) can be decomposed into the product AL^{kx)ALy{ky) with 
Aiik) — )■ 6{k) as L — )■ oo. Employing the decomposition of x given in Eq. (4.14) allows the 
density of the normal fraction to be related to the transverse component of x evaluated at zero: 

Pn = (3\unxt{k) = (3xt{0). (4.22) 

On the other hand [Fig. 4.1(c)] we may take the limit Ly oo first, resulting in an infinitely 
long channel — with velocity perpendicular to the walls — in which the entire body of the fiuid 
must move regardless of whether it is a superfiuid or not. In a similar way to the previous 
paragraph, pu = f3 limfc^_j.o ^^^ky^o x(k) ■ u, and making use of the decomposition in Eq. (4.14) 
we find that the total density is related to the longitudinal component of the correlations: 

p = (3]imxi{k) = ^Xi{0). (4.23) 
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With these expressions, the normal fraction may finally be expressed directly as 

fn = Pn/p = ]imxt{k)/]imxi{k) (4.24) 

while the superfiuid fraction is = 1 — fn- Thus, we have expressed the superfiuid and 
normal fractions in terms of a correlation function which can be directly computed in the 
thermal ensemble; there is no need to deal with difficult moving boundary conditions or other 
dynamical perturbations in the simulation itself. 



4.3 Numerical procedure 

To determine the superfiuid fraction from a PGPE simulation, we need to estimate the tensor 
of momentum density correlations x using simulation results. For a finite system constrained 
to a periodic simulation box — as studied in chapter 3 — we may only compute the momentum 
correlations at discrete grid points. The discrete analogue of Eq. (4.13) leads to the expression 

X(k) oc (pkP_k) (4.25) 

where the constant of proportionality is not important to the final result, and Pk are the discrete 
Fourier coefficients of p(x) over our simulation box. 
The momentum density operator is given by 



P(x) = y 



(V^'^(x))V'(x) - V'"f(x)V^/'(x) (4.26) 



which may be derived by considering the continuity equation for the number density, ('?/'^(x)'?/;(x)). 
For a given classical field Eq. (3.9), the Fourier coefficients of p may be written as 



= E(2k' + k)c*,Ck+k„ (4.27) 



where As is the area of the system. Computing a value for all pk at each time step, we then 
evaluate x(k) via the usual ergodic averaging procedure using Eq. (4.25). 

Having evaluated x(k), we are left with performing the decomposition into longitudinal 
and transverse parts. For this, simply note that Eq. (4.14) implies Xii^) = k ■ x(k) ■ k, and 
Xt{k) = w ■ x(k) ■ w, where w is a unit vector perpendicular to k. 

Values for xt and xi may be collected for all angles as a function of k, and a fitting procedure 
used to perform the extrapolation /c — )■ 0; this procedure is illustrated in Fig. 4.2. At low 
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temperatures, the extrapolation is quite reliable, but becomes more difficult near the superfluid 
transition where sampling noise increases and Xt{k) changes rapidly near k = 0. Without a 
known functional form, we settled for a quadratic weighted least squares fit of In(xt) and ln(x/) 
versus k. A weighting oi 1/k was used to account for the fact that the density of samples of x vs 
k scales proportionally with k due to the square grid on which is evaluated. The logarithm 
was used to improve the fits of Xt very near the transition where it varies non-quadratically 
near k = 0. The fitting procedure and extrapolation to A; = generally produces reasonable 
results, but is somewhat sensitive to numerical noise. For this reason, the computed superfluid 
fraction at high temperatures is not exactly zero (see Fig. 3.2 on page 49). 



1.1 - 




20 40 60 80 100 120 

k [dimensionless units] 



Figure 4.2: Example fitting and extrapolation to k = for the transverse and longitudinal 
components of the momentum density autocorrelation tensor, x- The apparent functional form 
for Xt and Xi changes with temperature — particularly near the transition — which along 
with the sampling noise makes them difficult to ffi reliably. The data shown corresponds to a 
temperature T ^ 0.99Tkt slightly below the transition. 

4.4 Discussion 

We have applied the technique described in this chapter to study the BKT phase in chapter 3. 
This included computing the superfluid fraction as a function of temperature — see Fig. 3.2 
on page 49 — and the shape of the curve is consistent with the expectation of a universal 
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jump^ in the superfluid density, as predicted by Nelson and Kosterlitz [92]. The location of 
the BKT transition — as deduced from the sharp disappearance of the superfluid fraction — 
is also consistent with the behaviour of the other physical quantities computed in chapter 3, 
including the decay of spatial correlation functions and the vortex unpairing. 

Unlike explicitly constructing a dynamical simulation, our method is based on correlations 
that are calculated in the stationary thermal ensemble. This is a great practical advantage 
because it avoids the need for additional numerically costly simulations. It also avoids any 
concern that the dynamical perturbation in such simulations might disrupt the thermal back- 
ground, resulting in a perturbed measurement of the superfluid fraction. Having said this, it 
would be interesting to compare our results to the superfluid fraction as deduced from the drag 
force felt by a moving impurity potential. 

We finish by discussing some limitations and possible extensions to our technique. The 
main limitation is that the method appears to rely critically on translational symmetry; at the 
very least it is clear that the thought experiment used in the derivation doesn't make sense for 
a trapped system. This is unfortunate, since all experimental systems are necessarily trapped 
in some way or another and any quantitative comparisons with actual experiments must take 
this inhomogeneity into account. So far the superfluid fraction for the trapped system has been 
determined by using the universality result for the critical density in the homogeneous gas [106] 
in combination with the local density approximation [68, 13]. It would be interesting to be able 
to compute the superfluid fraction independently as we have done here, but it is not obvious 
how to generalise the derivation to this case. 

In this thesis we have calculated the superfluid fraction only for the two-dimensional case. 
In ID the tensor is degenerate so there is no transverse component available, and the 
derivation fails to make sense. There is no such problem in 3D, and we expect our method to 
provide useful results in this case. 



^Note that we do not expect an exact discontinuity in our calculated superfluid fraction due to the finite size 
of the system and the presence of statistical noise. 



Chapter 5 

Effective ID equations for trapped 
Bose-Einstein condensates 



This chapter presents an ansatz for solutions to the 3D GPE in the quasi-lD regime. Our 
ansatz expresses the full 3D wavefunction in terms of a pair of ID complex fields that describe 
the amplitude and width of an elongated BEC. We derive equations of motion for the ID fields 
using the Lagrangian formalism and solve them numerically for several test scenarios. We 
compare with other ID approaches, and the full 3D solution. 

5.1 Introduction 

The Gross-Pitaevskii equation (GPE) described in section 2.3 [Eq. (2.68)] is a remarkably 
successful tool for describing experimental BEC dynamics near T = 0. The reasons are simple: 
it is accurate for a large subset of experiments, and efficient numerical schemes are fairly simple 
to implement. Nevertheless, solving the GPE in three dimensions can be computationally 
demanding because of the size of the spatial grid required. In some cases it is possible to 
mitigate this problem by dimensional reduction; there are at least two possibilities: First, we 
may make use of symmetries in the physical situation to reduce the number of dimensions. For 
example, experimental systems commonly have cylindrically symmetric trapping potentials; if 
the initial state also has cylindrical symmetry we may then eliminate the angle variable, which 
reduces the simulation to two dimensions. Second, there are experimental situations where 
the system is strongly confined in one or two "transverse" dimensions so that the dynamics 
in those dimensions is particularly simple. We may then make an ansatz for the wavefunction 
that allows us to integrate out the transverse directions to produce a lower dimensional effective 
equation which may be solved on a much smaller numerical grid. 

Elongated, cylindrically symmetric cigar-shaped BECs are routinely produced in the labora- 
tory using a variety of trapping techniques including optical lattices [56], atom chips [64, 79, 44, 
114] and other types of magnetic traps as first achieved by Gorlitz et al. [54]. In parallel, several 
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studies (see, for example, [34, 120, 65, 20, 85]) have investigated the dynamics of dispersive 
shock waves^ in BEC. Refs. [20, 85] combined both one dimensionahty and shock wave gener- 
ation by rapidly applying or removing a perturbing potential at the centre of a cigar-shaped 
BEC. In both of these studies the resulting shock waves were modelled numerically using a GPE 
or GPE-like mean-field approach. The calculations in Ref. [85] required a high spatial resolution 
in the longitudinal direction due to both the large size of the condensate and the need to resolve 
shock fronts at small length scales. As a result, these calculations were restricted to using a 
ID effective equation for efficiency. The method chosen was the non-polynomial Schrodinger 
equation (NPSE) of Salasnich et al. [Ill], one of several methods that adiabatically eliminate 
the transverse motion and assume the shape of the transverse profile varies slowly as a function 
of the longitudinal coordinate. 

While the NPSE is convenient, it is not obvious that the assumptions used in its derivation 
are valid when dealing with shocks. In particular, the existence of a shock implies rapid variation 
of the density and other system parameters with the longitudinal coordinate, and we would 
expect to this to carry over to the shape of the transverse profile. In this chapter we relax both 
of the assumptions which go into deriving the NPSE — we include both the variation in the 
longitudinal direction and avoid making the adiabatic approximation. 

An additional motivation for this work was to derive a ID effective equation capable of 
simulating the expansion dynamics after turning off the trapping potential. These type of 
expansions are the standard experimental tool for imaging condensates, but simulating them 
directly is difficult due to the large size of the spatial grid required. The alternative ID ansatze 
discussed in the next section eliminate the transverse velocity either implicitly or explicitly, and 
are therefore fundamentally unable to deal with expansion. 

5.1.1 Previous work 

We consider quasi-lD systems where the x and y coordinates correspond to the tightly con- 
fined transverse directions; the z coordinate is the weakly-trapped longitudinal direction. We 
assume that the transverse trapping potential is harmonic and cylindrically symmetric, while 
the longitudinal potential is generic. The full potential is 



dispersive shock wave occurs when dispersion rather than dissipation dominates the physics at small 
length scales. Dispersive shock waves are characterised by pulse trains as seen in the simulations of Ref. [85], 
for example. 




(5.1) 
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where r = ^ is the radial coordinate, x = {x,y,z), m is the atomic mass and uj± the 
angular frequency of the transverse harmonic potential. 

In the quasi-lD regime, uj± is large enough to prevent significant excitation of the transverse 
degrees of freedom. The simplest effective ID equation may be derived by considering the case 
without interactions, Uq = 0. In this case the full 3D wavefunction </> factorises and the 
transverse component has the functional form ^^7=^e~''^/^°"^, where a is the transverse width. In 
the perturbative regime^ where Uq is small but nonzero, this suggests the simple ansatz 

0(x,^)=V^l(^,^)^e-^'/2-^ (5.2) 

with a taken equal to the width of the non-interacting ground state. This ansatz works fairly 
well when Uq is small and leads to the ID GPE which describes the evolution of the unknown 
function tpi. 

thdt^i = -^dliJi + V^i + f/iDlV^iTV^i, (5.3) 

where Uid is the effective ID nonlinearity constant. For larger [/q the transverse width increases 
substantially and the ID GPE becomes rather inaccurate unless the value of a is adjusted 
accordingly. An appropriate value may be computed using a variational calculation, assuming 
a constant density in z. 

More sophisticated one-dimensional approximations have been investigated by several au- 
thors; we review those that are relevant to the current work below. An attempt has been made 
to keep to the notation used in the original papers, with some modifications for consistency. 

The first of the more sophisticated approximations was the NPSE, as described by Salasnich 
et al. [Ill] in 2002. They used the ansatz 

0(x, t) = ] e-V^-^(-'*)/(^, t), (5.4) 

along with the assumption that the width a changes sufficiently slowly in the longitudinal 
direction that dzcr is negligible. The equation for / which results is known as the nonpolynomial 
Schrodinger equation due to the nonpolynomial nonlinear term that arises after eliminating the 
transverse width a: 

^h^tf = -^^lf + V-J + 
2m 




■^The perturbative regime is defined by a^JT-i <C 1 where ni is the integrated ID density — see, for example, 
Ref. [89]. 
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where a± = ^/hjmu^ is the transverse length scale. Implicit in the ansatz from Eq. (5.4) is the 
assumption that transverse dynamics given by changes in a are much faster than the dynamics 
of / in which we are interested. To see this, note that the phase of does not depend on the 
coordinate r, which implies the absence of a radial superfluid velocity — the ansatz does not 
support transverse dynamics independently of the field /. 

This lack of transverse dynamics was addressed more carefully by Muiioz Mateo and Delgado 
[89, 90], who showed that the general ansatz 

0(x,t) = v9(r,ni(2,t))0i(2,t) (5.6) 

may be used to derive an equation for 0i via the adiabatic elimination of the transverse degrees 
of freedom. Here is some family of radial wavef unctions, parametrised by the local ID 
condensate density n\. By assuming that d^nx = 0, making the adiabatic approximation, and 
integrating out the transverse direction, they obtained the remarkably simple equation 

ihdt(j)i = -— (9^01 + V;0i + /i±(ni)0i, (5.7) 
2m 

where 

/i±(ni) = lldxdyip* ( -^Vi + + Uon,\ip\^\ ^ (5.8) 

is the local chemical potential. (We will refer to Eq. (5.7) as the Mufioz Mateo-Delgado equation 
(MDE) from now on.) Choosing the formula /i_L(ni) = /ioj^i/l + AagUi to interpolate between 
known limits at large and small Uq, they were able to obtain very accurate predictions of in- 
trap oscillations with a range of smooth initial conditions, when benchmarked against a full 3D 
calculation [89]. 

When the transverse dynamics are important or the condensate width varies rapidly with 
the longitudinal coordinate 2, we expect approaches based on the adiabatic approximation 
to become less accurate. Kamchatnov et al. [73] investigated a different generalisation of the 
NPSE ansatz, 

Vvr6(2;,t) 

using a variational formalism. The added generality and lack of assumptions regarding the z 
derivatives of the fields were motivated by the desire to describe solitons involving short length 
scales. The Kamchatnov ansatz results in coupled equations^ for the two real fields 6(2;, t), 0(2;, t) 
and one complex field 'ijj{z,t), which were used to analytically investigate solitonic solutions in 



■^We refer the reader to Ref. [73] for details of the equations — they are not used further here, and are 
reasonably complex to write down. 
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various limits, in addition to small amplitude linear waves. An important feature is the inclusion 
of the phase factor a, which allows for "breathing mode" dynamics in the transverse direction. 

In this chapter we introduce a variation of the Kamchatnov approach, which — while 
mathematically equivalent — leads to more compact equations. In contrast to the earlier work, 
we investigate general numerical solutions. We show cases where the solutions are accurate, and 
highlight a number of generic difficulties which will arise for any similar factorisation ansatz. 

5.2 Formalism 

5.2.1 The Ansatz 

We consider the following approximation to the 3D wavefunction: 

0(x,t) =^(^,t)e^("'*)"'. (5.10) 

Here both ip and x complex time-varying fields in one spatial dimension. The field x is 
related to the transverse Gaussian width and transverse superfiuid speed, respectively, via 

a{z) = 1/a/— 2Rex(^) and v±{z,r) = 2\lmx{z)\r. (5-11) 

The field tp contains density and longitudinal phase information along with a normalisation 
factor for x; with this ansatz, the ID density has the form 

„(M)^//^.<<.Wx,*)P^3w^fffj;;^. (5.12) 

In what follows we will often suppress time and space arguments to avoid excessive notational 
clutter. 

5.2.2 Equations of motion 

To derive the equations of motion we use a time dependent variational formalism essentially 
the same as that described in section 2.3.1. We begin by computing a Lagrangian for the one 
dimensional effective theory by integrating out the transverse dimensions. Taking derivatives 
of the effective Lagrangian and putting these into the Euler-Lagrange equations then yields the 
equations of motion in the usual manner. 
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One possible Lagrangian^ for the three-dimensional classical |0| field theory is 



ih 



(5.13) 



where the Hamiltonian functional is given by 



(hi. 



2m 



(5.14) 



Differentiating C with respect to </>* yields the 3D GPE via the Euler-Lagrange equations. If 
instead we substitute our ansatz from Eq. (5.10) into Eq. (5.13) we may perform the Gaussian 
integrals over the transverse coordinates, x and y. We then obtain a reduced Lagrangian which 
is a function of the fields ip and x'- 



dz 



ih 



x-x* 



2 lix + X*)^ 



x + x* 



2m [{x + X*y 



ix + x*)^ x + x* 



+ 



^d^ip^d^X + iJ*d^i)d^x* 



+ 



mu 



{x + x*Y 

2 L/,|2 



X + X* 2(x + X*) 4(x + X*) 



(5.15) 



As usual, stationarity of the action implies the Euler-Lagrange equations for our two fields, 
which are 



d 5C 



5C 



d 6C 



5C 



dt\6x*J Sx* dtyStp*) {^-"^Q) 

Performing the functional derivatives (see appendix A), rearranging and simplifying the results 
gives the equations of motion for x ^iid ip: 



2m 



x + x* 

+ yA + 



+ 



muj\ 



2 



+ ^I^P(X + X*), 



(5.17) 



We highlight two relevant features of the equations above that may not be immediately obvi- 
ous: First, Eqs. (5.17) are energy conserving due to the time independence of the Lagrangian^. 
Second, our choice of generalised coordinates in Eq. (5.10) forces us to give up the Hamiltonian 
structure of the phase space in exchange for a compact parametrisation of the wavefunction. 



^This form of the Lagrangian is manifestly symmetrical with respect to and 0* but we could equally well 
have chosen a form more similar to that used in section 2.3.1. 

^The sophisticated way to say this is that the time translation symmetry of the Lagrangian implies energy 
conservation via Noether's theorem. 



5.2 Formalism 



77 



This is not physically problematic, but does affect our choice of numerical methods. 



5.2.3 Conservation of normalisation and conserved current 

An important property of any low- dimensional effective equation arising from the GPE is the 
conservation of normalisation of the wavefunction; this corresponds to conservation of the total 
number of atoms during time evolution. The total normalisation is 



N 



(hi. 



dz 



-(x + x*)' 



(5.18) 



In principle, is a function of time; to prove that it is not and the total number is 
conserved, we show that = 0. Expanding A^ and inserting the equations of motion to remove 
the resulting factors of if) and we have: 



N = -Ti I dz 

ihn 



2m 



dz 



x + x* (x + x*)^ 

(x + x*)^ 



X + X* (x + x*)^ (x + x*)^ 



— c.c. 



(5.19) 
(5.20) 



where c.c. stands for the complex conjugate of the preceding bracketed term. Finally, we may 
remove all second order derivatives using integration by parts, assuming that the boundary 
terms are zero^. The resulting terms cancel out, implying that A^ = and normalisation is 
conserved. 

Conservation of normalisation suggests that there should also be a conserved one-dimensional 
current j{z,t), obeying the equation 



h{z,t) = -dj{z,t). 



(5.21) 



Computing n directly from the definition in Eq. (5.12) and rearranging shows that the appro- 
priate conserved current may be written 



3 



h 

2mi 



(x + x*)^ x + x* 



— c.c. 



(5.22) 



If we define the ID wavefunction r] = 4'\/—tt/{x + X*) with the properties n{z,t) = |?7(z,t)p 



®In a periodic system the individual boundary terms are generally nonzero but nevertheless correctly cancel 
due to periodicity. 
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and arg [0(0, 0, 2;, t)] = arg [77(2, t)], the current takes on a somewhat more famihar form 

J = 7^ (v*d.v - Vd.v* + \V\'^^^^^-^) ■ (5.23) 
2mi V X + X J 

The usual conserved probabihty current for quantum mechanics consists of the first two terms 
in the expression for j. The extra term is a consequence of integrating out the radial structure. 



5.2.4 Ground states 



We make use of ground states as physically reasonable initial conditions for the numerical 
simulations in the next section. Here we describe the ground state equations to be solved, 
along with the simplest approximate solution. 

Ground states of Eqs. (5.17) have no transverse dynamics, that is, x = 0- Further, the 
transverse speed v± is zero which implies Imx = 0. The time evolution of tp is the simple 
phase rotation ip = {—i^/h)il) for some chemical potential and we may assume for simplicity 
that ip is real at t = 0. With these observations, we see that the ground state obeys the 
time-independent system of equations 







1_ 

2m 



0^^ 

2m 



-dlx 



2M5.X + 2^-4x^ 



{dzxYi^ 

2x' 



X 



- ^i'x 



^ ^ /|2 

3t/o, 



(5.24) 



where both x and ip are real. 

In cases where the local terms dominate in Eqs. (5.24) — slow spatial variation, large 
interactions, or large densities — we may ignore the z derivative terms. This is equivalent to 
using X and ip as determined from a spatially homogeneous system with the same local density, 
and is therefore known as a local density approximation (LDA). 

Setting dzip = and d^x = 0, Eqs. (5.24) reduce to 



2hl_ 

3m 

„;.|2 



x' + 7T(/i-"^)X + 



A 



0, 



3t/n 



2h^ 



m 



X + (/i-^) 



(5.25) 



These algebraic equations are easily solved for x and ip. We take the branch of the square root 
such that X < so that the transverse width in Eq. (5.11) is a real number. We note that 
Eqs. (5.25) are analogous to the well known Thomas-Fermi approximation for the GPE ground 
state. 
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5.3 Numerical Simulations 

In this section we present numerical results for the simulation of Eqs. (5.17) along with com- 
parisons to full 3D simulations and other proposed effective ID equations. Our ID equations 
are discretised on a uniform grid with periodic boundary conditions; derivatives are calculated 
spectrally via the Fourier transform. We use the standard fourth order Runge-Kutta integrator 
for time stepping. For the numerical simulations we use units such that h = m = 1, and have 
arbitrarily chosen u± = 10 and a periodic domain —8 < 2; < 8 with wavefunction normalisation 
= 1000. With these choices, nonlinear behaviour becomes strongly apparent at values of the 
interaction strength Uq ^ 1. 

The full 3D system was assumed to be cylindrically symmetric, allowing a reduction to 
a two-dimensional equation. We simulated this 2D equation using the XMDS2 software [1] 
with a Bessel Fourier basis. Time stepping was achieved using an adaptive fourth/fifth order 
Runge-Kutta solver^ for time stepping. 

5.3.1 Ground states 

A common approach for finding the ground state of the GPE is to minimise the Hamiltonian 
via continuous-time steepest-descent optimisation, the so-called "imaginary time" method*. 
Unfortunately, this method does not always work well for the equations presented here, because 
the energy depends only weakly on the transverse degrees of freedom x in regions of low density. 
This means that convergence of x to the ground state value can be extremely slow in regions 
of high potential V^. 

To avoid this problem, we solve a discretised version of the ground state differential equations 
(5.24) directly using Newton's method. The spatial discretisation and derivative calculation 
used here is the same as for the dynamical equations. The desired normalisation for the wave- 
function is obtained by treating the chemical potential as one of the unknowns and adjoining 
the normalisation condition to the set of equations. 

Convergence of Newton's method for this system requires a starting point which is quite 
close to the true solution. For this we make use of the LDA solution Eq. (5.25). We find that 
the presence of low density regions where the LDA gives ip = results in Newton's method 
failing with a singular Jacobian. A simple method to avoid this problem is to smooth the LDA 
solution with a spatial filter with exponentially decaying tails. 

''named ARK45 in the XMDS software 

*We note that for the equations presented here, steepest-descent minimisation is not equivalent to replacing 
the time t with ir and evolving in r. This is because the equations of motion are not Hamiltonian after the 
transformation in Eq. (5.10). 
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We also needed to find ground states for tlie NPSE, MDE and 3D equations for comparison 
purposes. The imaginary time method is sufficient in these cases because the equations arise 
directly from 3D or effective ID Hamiltonians, albeit with sometimes unusual nonlinearity 
terms. To satisfy a given normalisation we use a slightly unusual version of the imaginary time 
method where the chemical potential is treated as an unknown and adjusted continuously 
along with the wavefunction (see appendix C for details). 

5.3.2 Test case: soliton formation 

As a test case, we consider evolution in a circular waveguide {Vz = 0) after releasing the ground 
state of the barrier potential Vz = C exp(— lO^^) at time t = 0. The barrier intensity C = 50 is 
chosen such that the density inside the barrier is depressed to few percent of the background 
density to induce strongly nonlinear evolution. 

After finding the ground state, the barrier is turned off and the atoms fill the resulting hole. 
For sufficiently large interaction strength Uq, the excitations are carried away as a combination 
of sound waves and a pair of grey solitons form, as shown in Fig. 5.1(a). The associated phase 
profile presented in Fig. 5.1(c) shows the expected jump in phase for a soliton across the density 
depression at t = 1. In addition there is a decrease in the width and associated radial flow. 

To evaluate the ID model we compare with the full 3D results, and also to two alternative 
effective ID equations: the NPSE and the MDE. A comparison of the soliton speed as a function 
of interaction strength is shown in Fig. 5.2(b). We see that our equations and the MDE are 
competitive over most of the range, however the MDE wins out at larger interaction strengths. 
Of note is the NPSE results that are significantly less accurate than our equations, even though 
the same transverse Gaussian ansatz is used. 

We note that the MDE has a particular advantage in the reproduction of accurate ground 
states: the interpolating form used for the local chemical potential accurately represents the 
important properties of the transverse wavefunction. In contrast, both methods which assume 
a transverse Gaussian profile are less accurate when the transverse wavefunction is deformed 
due to the interaction energy. Nevertheless, our equations are significantly more accurate than 
the NPSE, which indicates that derivatives of the transverse wavefunction parameters have a 
measurable effect for these initial conditions. 

Further examination of Fig. 5.1(a) shows an apparent oscillation in the soliton depth coupled 
with the emission of sound waves. This feature is not present in the MDE or NPSE simulations 
but is present in the full 3D results. To quantify the effect. Fig. 5.2(a) shows the soliton depth 
as a function of time. We see that the initial condition causes an oscillation in the transverse 
width which feeds back into the soliton depth in this case. While the width cannot be observed 




Figure 5.1: Solitons on a constant density background. The initial condition is the ground 
state of a Gaussian barrier potential, (a) Spatial density as a function of time, (b) Density 
slice through the centre of the inferred 3D cloud at t = 1. (c) Detail of the fields at t = 1 — 
dashes indicate the phase profile, aig^ip), the solid line is the width of the Gaussian, a, and the 
dotted line shows the scaling factor for the radial velocity, Im x- 
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Figure 5.2: Comparison between various simulation methods for the sohton waveguide problem 
of section 5.3.2. (a) The soliton depth as a function of time for Uq = 0.25. The high frequency 
oscillations in the 3D solution correspond to a small excitation of the higher order transverse 
modes, (b) The soliton speed as a function of the interaction strength Uq. 
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after integrating out the transverse direction, the sohton depth can. 
5.3.3 Transverse shock wave formation 

We have presented a case above where our equations are competitive with other effective ID 
equations for modelhng the full 3D dynamics. Nevertheless, we have found recurring stability 
problems over the course of many numerical experiments with varying initial conditions. Such 
instability is typically characterised by sharp shock-like features that initially develop in the 
transverse width x- We observed that the shocks are often associated with the density passing 
close to zero, and inspecting the equation of motion for x Eqs. (5.17) suggests that the 
term containing l/ip might be responsible: When ip is very small this term becomes very large, 
potentially resulting in stiff equations. Explicit schemes such as the classic fourth order Runge- 
Kutta method become unstable in the presence of stiffness unless an unmanageably small time 
step is used. 

We tested this hypothesis by making use of a Matlab implementation [38] of the RADAU5 
algorithm [63]. RADAU5 is a fifth order instance of the Radau IIA class of implicit Runge- 
Kutta algorithms, and is known for its exceptional stability properties [63]. Nevertheless we 
found the results to be inconclusive: while it seemed to help for some initial conditions, there 
were certainly cases where the solutions diverged regardless of the type of numerical integrator 
used. 

After further investigation, we believe we have identified the underlying cause of these 
problems: The ideal transverse width is not always a continuous function of the transverse 
slices of the 3D wavefunction. To make matters worse, discontinuities develop dynamically and 
it is not generally obvious which initial conditions will lead to problems. 

To demonstrate the issue, we fix Uq = 0.1 and consider initial conditions that are the ground 
state of the trapping potential 



By changing Vq, the nonuniformity of the initial conditions can be increased from a completely 
homogeneous system (Vq = 0) to one in which all the density is concentrated near z = 0. For 
sufficiently strong initial nonuniformity (Vq = 50) we observe the formation of discontinuities 
in the transverse width of the 3D system at a finite time t ^ 0.17 after the pulse is released, as 
shown in Fig. 5.3. We computed the 3D width by nonlinear least squares fitting of a Gaussian to 
the 3D transverse density profile shown in Fig. 5.4. Also shown in Fig. 5.3 is the width arising 
from our ansatz, according to Eq. (5.11). While the width oscillations are not reproduced 
accurately, the existence and position of the discontinuity is correct. As a result, numerical 
methods that assume x is smooth will generally fail within the next few time steps after the 



V{z) = Vo 1-e 



(5.26) 
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snapshot shown in Fig. 5.3. 




Figure 5.3: Transverse Gaussian width (standard deviation a) of the wavefunction density from 
our ansatz, compared to a 3D simulation at time t = 0.17. This is a zoomed view in the z axis, 
clearly showing the discontinuities at 2 ~ 1.6. 



The existence of discontinuities comes as a surprise; one might expect that continuous 
evolution of a smooth underlying field should lead to smooth fitting parameters. The 3D 
transverse profile in Fig. 5.4 shows that the underlying density is smooth as expected, and that 
the transverse profile is well behaved along most of the z coordinate. However, in the region of 
very low density near the discontinuity {z ^ 1.6), the transverse profile has two peaks as shown 
in Fig. 5.5. The discontinuity comes about because the least squares solution jumps suddenly 
from fitting the thin central peak to fitting a wider combination of the central and secondary 
peaks. Even though the density is low in the problematic region, the discontinuity in width is 
enough to destabilise the numerical solution. 

We expect that any ansatz for the 3D wavefunction containing a single transverse width 
parameter will suffer from similar problems. It is somewhat ironic that the NPSE, MDE and 
similar approaches based on the adiabatic approximation remain stable and somewhat accurate 
even in regimes where our equations fail: Such methods avoid the problem by assuming that 
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z 

Figure 5.4: Density snapshot from a 3D simulation showing radial structure at time t = 0.17, 
as in Fig. 5.3. 

the transverse profile depends only weakly on z\ 

5.4 Conclusion 

In this chapter we have examined the general Gaussian ansatz for the quasi-lD Bose gas tightly 
confined in a transverse 2D harmonic trap. We have compared the ansatz to two alternative 
ID effective equations, the NPSE and MDE. We have found that some transient behaviours — 
such as an oscillation in soliton depth after release from a trap — can only be captured with 
an ansatz that allows transverse breathing motion. Such motion is prohibited by the NPSE 
and MDE due to an assumption that the transverse degrees of freedom adiabatically follow the 
local density. Further, we have shown that the method is more accurate than the NPSE due 
to the inclusion of derivatives of the transverse width. 

However, we have found that our general Gaussian ansatz has several important failings 
that make it impractical for general use. First and most severely, discontinuities appear in the 
Gaussian width parameter for a wide range of sufficiently excited initial conditions. Second, 
finding ground states numerically is significantly more difficult than with competing ID ap- 
proximations. Third, the assumption of a Gaussian profile implies less accurate ground states 
than the MDE when the interaction strength is sufficiently strong. 

The appearance of discontinuities in the width at finite time is an interesting issue which 
we trace to the appearance of two peaks in the transverse structure. Although this problem 
typically occurs in regions of low density and energy, it spells disaster for numerical solutions. 
We expect it is a generic feature of any similar ansatz which fits the transverse wavefunction 
using a single width parameter. One could consider additional parameters to get a better fit 
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Figure 5.5: Two transverse Gaussian fits (red) to the radial profile of the 3D wavefunction 
(blue) at adjacent grid points in the z direction. This example shows how a discontinuity in 
the fitting parameters may arise even when the underlying function to be fitted is continuous. 
The grid points in z span the position of the discontinuity in a shown in Fig. 5.3. 

to the transverse wavefunction, such as taking the next one or two harmonic oscillator basis 
states, scaled by a Gaussian width parameter. However, the resulting equations would clearly 
be more complex to derive — and more unwieldy to use — than those presented here. 

Our results show that out of the three effective ID equations considered, the MDE [89, 90] 
has the best combination of reliability, simplicity and relative accuracy. 




Chapter 6 

Conclusion 



In this thesis we have used numerical simulations to investigate the physics of ultracold Bose 
gases confined to one- and two-dimensional geometries. 

Our two-dimensional simulations were carried out using the PGPE classical field technique. 
The PGPE allowed us to investigate the finite temperature regime surrounding the BKT su- 
perfiuid phase transition, in a parameter regime that is relevant to current experiments. We 
directly computed the superfiuid and condensate fractions and found that these drop to zero at 
the same temperature: in our finite-sized system there is no clear separation between BEG and 
BKT transitions. This is in contrast to what is expected based on the theory of the infinite-sized 
system, emphasising the importance of finite-size effects. 

The unbinding of vortex pairs above the BKT transition is one of the defining features of the 
transition, and we developed a coarse-graining technique to quantify the number of unpaired 
vortices in the system. Results proved consistent with several other characteristic features, 
including the decay of the first-order correlation function and the sharp drop of the superfiuid 
fraction to zero. 

With a detailed microscopic theory at hand, we were able to simulate and validate the 
experimental measurement technique of Ref. [60], including the method used to deduce the 
first-order correlations from interference experiments. In addition, our simulations suggest that 
direct observation of vortex pairs is not possible due to limited imaging resolution. However, for 
this same reason it is likely that the dislocations observed in experimental interference patterns 
correspond to free vortices, and accordingly are a strong indicator of the BKT transition. 

To investigate the BKT transition we required a method for calculating the superfiuid 
fraction from a PGPE classical field simulation. Ideally such a method would be based on 
properties of the equilibrium thermal ensemble as ergodically sampled by the PGPE. We were 
able to derive such a method, using linear response theory to relate the superfiuid fraction 
to the autocorrelations of the momentum density operator. Our method provided physically 
reasonable results that were consistent with the other physical properties of the 2D gas. 

One-dimensional effective equations allow for efficient simulation of elongated cigar-shaped 
BEGs. We derived equations of motion using the Lagrangian formalism with a Gaussian ansatz, 
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and solved these numerically for several test cases. We found that in certain cases our equations 
are more accurate than effective equations based on the adiabatic approximation. However, 
they are numerically unstable for a wide range of initial conditions. We investigated this 
instability and found that it is an inherent weakness of any similar ansatz which relies on a 
single variational width parameter. 

In summary, this thesis contributes original insight into the physics of ultracold Bose gases 
in low dimensions, shedding light on recent experiments and supplying several new theoretical 
and numerical tools. 
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Appendix A 

Mathematical techniques for classical field 

theory 



A.l The Wirtinger calculus: derivatives of nonholomor- 
phic functions 

The Wirtinger derivative is the suitable generahsation of the real derivative to general complex 
valued functions for the purposes of making linear approximations. As we will see below, the 
usual derivative of complex analysis is much too strict for these uses. 

The derivative in complex analysis is defined using the familiar limit definition used in the 
real calculus, only with complex variables. Consider a function /: C — C and some point 
zo E C; the derivative of / at zq is 



(see, for example, [22, page 30]). If the limit converges to the same value irrespective of the 
way in which z zq we say that / is differentiable at ^o- A function is called holomorphic on 
some domain D if it is differentiable at every point in an open set containing D. The existence 
conditions for complex derivatives defined like this are much stricter than for the real derivative. 
In particular, a necessary condition for the existence of the limit in Eq. (A.l) is the Cauchy 
Riemann conditions: If f{x + iy) = u{x, y) + iv{x, y) with u, f : — )• R, the Cauchy- Riemann 



Deriving the equations of motion for a classical field theory from an action principle depends 
on taking derivatives of the action functional with respect to the fields and setting these to zero. 
That is, for an action S[%1)\ we would like to do something like set 




Z- Zq 



(A.l) 



conditions are = and |^ = -fi. 

ox ay ' ax ay 




(A.2) 
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which should give us the equations of motion. However, the action S is generally a real function 
for any complex value of the field ip, and as such cannot possibly satisfy the Cauchy-Riemann 
conditions^ . From the point of view of traditional complex analysis this immediately leads to a 
problem: how do we make use of the derivative 6S/6ip when this derivative isn't even defined 
according to (some functional generalisation of) the definition in Eq. (A.l)? This is the main 
problem that we tackle in this section. 

One complicating factor for our particular case is that the derivatives of interest are func- 
tional derivatives with respect to the field ip. Luckily, the central problem remains intact if we 
forget about the fields and consider a nonholomorphic function / of a single complex variable 
z. This simplification is made throughout the rest of the section. 

If we can't use the usual complex derivative for common problems, what is the alternative? 
Derivatives are all about linear approximations^, so to understand what to do about the peculiar 
strictness of the Cauchy-Riemann conditions we consider a linear expansion of / about the point 
^0 = Xo + iyo in terms of the real variables x and y. Assuming that u and v are differentiable, 
this is simply 



f{x + iy) 




dv 
dy 



Ay 



(A.3) 
(A.4) 



where Ax = x — xq and Ay = y — y^. This is a straightforward expansion using the real 
calculus so it's clear that this linear approximation is valid regardless of the Cauchy Riemann 
conditions. Using Az = z — zq, Az* = z* — Zq for the conjugate and suppressing the subscripts 
for brevity (all derivatives being evaluated at zq), we have 



f{z)^f{zo) + ^Ax+^Ay 

^/(.o) + g^(A. + A.*) + |l(A.-A.*) 



dl 
dx 



dy. 



dl 
dx 



dy 



Az* 



(A.5) 
(A.6) 
(A.7) 



This is an elementary but interesting result: the general linear approximation to any real- 
differentiable complex-valued function may be expressed as a linear combination of Az and 



^At least not in a nontrivial way: for a real function / we have v{x, y) = and the derivative conditions 
immediately show that u{x, y) must be a constant. 

■^In fact, the Frechet derivative — an elegant generalisation of the derivative to Banach spaces — is by 
definition exactly the linear part of the affine approximation to a function. 
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A2;* with the coefficients given by a simple combination of the real derivatives. Furthermore, 
we can now express the ratio inside the limit in Eq. (A.l) as 

/(.)-/(.„)_! /a/ 1/9/^,9/ V-.„ (^3) 



A2; 2 \dx dy J 2 \dx dy 

where 6 = arg A2; which reveals the general form of the directional dependence of the limit. If 
/ is to be complex differentiable, the value of the right hand side must be independent of z 
which implies ff + ^f^ = 0; rearranging we see that this is exactly an expression of the Cauchy 
Riemann conditions. 

Coming back to Eq. (A. 7), we define two new derivative operators 

df ^1 /df ^df\ df _1 /df .df\ 

dz 2 \dx dy J dz* 2 \dx dy J ' 

which allow us to write the linear approximation as 

/W-/(^o) + |fA. + ^A.*. (A.10) 

At first sight this expression may seem like a mere curiosity that has been manufactured to 
resemble the real calculus of functions of two variables. On the contrary, it turns out that 
there is a powerful and general calculus for the operators ^ and Most importantly, given 
a function f{z) written in terms of z and z*, it is possible to compute derivatives directly 
with the new operators using the familiar rules of calculus for two real variables, and without 
decomposing / into real and imaginary parts. This calculus is known as the Wirtinger calculus 
[109] after the Austrian mathematician Wilhelm Wirtinger^. 

A. 1.1 Properties of the Wirtinger derivative 

We now discuss the general properties of the Wirtinger derivative as discussed in Ref. [109, 
§1.4]. Several properties are immediate from the definitions in Eq. (A.9): 

• The derivative operators are linear, due to the linearity of the real partial derivatives. 

• They obey the usual product rule. 

• The behaviour under complex conjugation is given by 



^It has also been called the CR-calculus by Krcutz-Dclgado [76]. 
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The derivatives satisfy the "independence of z and 2;*" property, 



dz* ^ dz . . X 

— = and 7^ = 0. A.12 

oz oz* 



The last property has great importance from a practical computational viewpoint, and we 
return to it after discussing the chain rule. 

Using the conjugation property and the linear approximation formula in Eq. (A. 10) it is 
straightforward to show that the Wirtinger derivatives obey a familiar looking chain rule. We 
have 

f{9{z)) ^ f (gizo) + ^Az + ^AzA (A.13) 



ft ( ,'9fdg df dg*\ . ^ f df dg df dg* 



dz dz 
d£dg_ ^ I 

dg dz dg* dz J \dg dz* ' dg* dz* 

where |^ = f^" have used the linear approximation formula twice, first for g and 

then for /. Identifying the coefficient in front of the Az using the linear approximation formula, 
we have shown that the chain rule for Wirtinger derivatives is 

^ ff ( w . df dg* 

d-z^^'^'^^=d-gd-z^W^' ^^-''^ 

with the analogous expression for -^f(g(z)). It is worth noting the close resemblance of this 
chain rule with the formula from the more familiar calculus on R^, 

^f(g(x, y)) = + (A.16) 

dx ' dgi dx dg2 dx 

where f,g:R'^^R^ and g{x,y) = [gi{x,y), g2ix,y)]. 

The relationship with the complex derivative in Eq. (A.l) of a holomorphic function is 
simply that 

the definitions coincide for this important special expected from the linear approxi- 

mation formula. For the same reason, it's possible to see that the Cauchy-Riemann conditions 
have a particularly nice form in this formalism: 
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The generalisation of the Wirtinger calculus to higher dimensions is completely straight- 
forward from a calculational viewpoint, because the variables Zk are genuinely independent 
for distinct k. We summarise this here by asserting that everything works just as in the real 
multivariable calculus, with the main difference being that variables always occur in conjugate 
pairs. 



A. 1.2 Computing with Wirtinger derivatives 

The independence property ^ = — along with the other basic properties — leads to the 
commonly repeated^ mnemonic device 

To compute with the Wirtinger derivatives, treat z and z* as "independent variables" 
and proceed exactly as for the calculus of two real variables. 

This "independence" is clearly problematic in a purely algebraic sense, but is extremely use- 
ful for the purposes of computing derivatives because it results in all the right rules when 
interpreted correctly. It is often useful to write functions f{z) redundantly as f{z,z*) to aid 
in applying the mnemonic. For example, f{g{z)) would become f (^g{z, z*), g*{z, z*)) and the 
chain rule follows after remembering the version for two independent real variables. 

We now present two examples that illustrate how computing with the Wirtinger calculus 
avoids the need to break a function into real and imaginary parts. A simple but useful example 
is the derivative of the absolute value: 

where the second line is an application of the chain rule with the square root as the outer 
function and noting that the second term that appears in the chain rule vanishes because 
-^z^^"^ = 0. In a similar way, = 

As a more complicated example, suppose we wanted to find the stationary points of the real 
valued function 

/(^) = |(^ + 1)^° + ^*|. (A.22) 

(We present this example because it is closely related to finding the stationary points of an ac- 
tion which is generally also a real valued function of complex arguments.) For general complex- 



*See, for example, Ref. [53, page 598] 
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valued functions, the stationary points are the points at which ff = and ^ = simulta- 
neously, but for real valued functions we have ^ = and it suffices to satisfy the first 
condition. Computing the derivative is an exercise in applying the chain rule to h{g{z)) with 
functions h{z) = \z\ and g{z) = {z + + z*: 



df dh dg dh dg* 
dz dg dz dg* dz 



(A.23) 



2\{z + iy^ + z*\J ^ ' \2|(2+l)io + ^ 

10 {{z* + 1)10 + z] (2 + Xf + (2 + 1)10 + z 




(A.25) 



2|(z + 1)10 + 2*1 

Therefore, finding the stationary points in this case involves solving the high order polynomial 

10 [(2* + 1)10 + 2] (2 + 1)0 + (2 + 1)10 + 2* = 0. (A.26) 

This could also have been derived by finding and setting hoth real derivatives to zero, but doing 
so using the formalism of the Wirtinger calculus is more convenient in this case. 



A. 2 Functional derivatives 

In the Lagrangian formulation of dynamics the action S* is a functional — that is, a function 
taking other functions as arguments^. Stationarity of the action means that the derivatives 
with respect to the fields are zero; this section investigates what we mean by such derivatives 
and how to calculate them. 

Throughout the section we make two simplifying assumptions to streamline the presentation: 
First, the functionals under consideration and their arguments are real. Second, the argument 
fields are functions of a single independent variable so that all integrals are one-dimensional. 
Lifting these assumptions presents no particular difficulty and is covered very briefiy in the 
next section. 



A. 2.1 An example 

As in the previous section, we start by looking at linear approximations. Consider for example 
the simple functional 

S[f]^ J dxg{x)f\x), (A.27) 

^This is obviously true for field theories, but remains true when we have a finite number of dynamical 
variables because such variables are functions of time. 
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where g is some fixed function, and / may be anytliing tliat is well behaved enough so that the 
integral is finite. Evaluating S at f{x) + t][x) for some small test function rj and expanding 
around /, we have 

S[f + V] = J dxg{x)[f{x) + vix)Y (A.28) 
dx g{x) [f{x) + 2f{x)ri{x) + ri^{x)] (A.29) 
dxg{x)f{x)+ dx 2g{x)f{x)r]{x), (A.30) 



where we have discarded the quadratic terms in rj. We see that the change in S which results 
from changing / to / + r/ is approximately a linear functional: 



S[f + V]- S[f] ^Lf[7]]^ J dx 2g{x)f{x)v{x). (A.31) 

In some sense the linear functional L defined here actually is the derivative of 5 — known as 
the Frechet derivative — but a slightly different definition is usually favoured by physicists, 
based on an analogy with the gradient. 

To motivate the definition of functional derivative commonly used in physics we compare to 
the discrete setting. For a finite number of variables the analogue of S would be some function 
i?(f) = J2iLi9ifi ^'^^ f ! g ^ R-^- The linear approximation to R about f is 

R{f + e) - R{f) ^ Yl ^9ifiei = Yl = Vi? ■ e. (A.32) 

i=i i=i •^^ 

The correspondences with the continuous case are / of, ?7-H-e, (7f-)-g and x ■H- i. With 
this in mind, it's clear that the function 2g ■ f in Eq. (A.31) is acting in an analogous role to 
the gradient V-R = [§f', • • • ? §f~] — [2fi'i/i; • • • ? '^Qn/n]- The functional derivative is therefore 
written using similar notation to the partial derivative: 



(A.33) 
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A. 2. 2 General definition 

Roughly speaking, the functional derivative of a functional 5* evaluated at / is a distribution^ 
4f such that 

S[f + S[f] = jdx ^^Vix) + O(hlp) (A.34) 

for all test functions rj. (Note that a more formal and general definition may be given in 
terms of the Frechet derivative^.) As indicated in the previous section, the distribution |j 
is a generalisation of the gradient vector from a finite to an infinite number of independent 
variables. In the case that the distribution is simply a function, we may evaluate it at x to give 
the number -^j^ which is intuitively the answer to the question "how much does S change if 
we change / by a small amount at the point x?" . 

This intuition leads to the "physicist's definition" of functional derivative: 

5S ^ S[f + e6..] - S[f] 
5f{x) e^o e 

where S^iy) = S{x — y). The notion of poking / exactly at x using a delta spike is intu- 
itive but this definition doesn't make much mathematical — or even calculational — sense if 
taken literally. For example, consider attempting to calculate the derivative of the functional 
J dx g{x)f'^{x) from the previous section. Trying to compute >S'[/ + e6x] then results in a term 
containing the integral of 6^ which has no well-defined value. It is possible to repair these 
problems by considering a sequence {A^ ,i} of nascent delta functions such that Ax^n — ^ as 
n — oo, and taking the limit e — )■ before the limit n — > oo. However, it seems better to simply 
regard Eq. (A. 35) as an aid to the intuition, and fall back to the idea of linear approximations 
for computational purposes. 



A. 2. 3 Calculating with functional derivatives 



Functional derivatives may usually be calculated directly by considering the linear expansion 
of the functional, but in practice this can be tedious. Therefore it is useful to derive a few 
general formulae that can be applied to a wide range of functionals. Here we present a few such 
formulae that are relevant to this thesis, along with an indication of how to derive them. We 



^ A distribution is also known as a generalised function and includes such objects as the Dirac delta "function" ; 
see, for example, Ref. [69]. 

''Borrowing the definition of the Frechet derivative, we might define |j as the distribution that satisfies 



hm 

MHO 



S[.f + S[f] - I dx 



ss 

Sf{x) 



r]{x) 
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direct the reader to Ref. [57, §2.3] for additional exposition. 



The functionals of interest are generally of the form 



S[f]= / dy g{f{y)J'iy)) 



(A.36) 



for some differentiable function g of two variables. The task is to create a linear approximation 
for 

S[f + r]] = J dyg{f{y) + r^{y)J'{y) + v'{y)). (A.37) 

Writing g^^''^^ and g^'^'^^ for the first derivatives with respect to the first and second arguments 
respectively, and briefly dropping the explicit y dependence of all functions for succinctness, 
the linear approximation is 



S[f + V] 



dy 



S[f]+ I dyg^'''\fj')r^+ I dy g^''^\fj')^' 



JO,l)(f 



(A.38) 
(A.39) 



This is almost in the right form, but the last term has the derivative t]' of the test function 
rather than the test function itself. The derivative may be removed by integrating by parts 
and assuming that the boundary terms are zero. This is true for a periodic domain, an infinite 
domain where the fields and their derivatives decay at infinity, or any domain where the test 
function t] is constrained to be zero at the boundary. Integrating by parts leads to 



S[f + v]^S[f] 



and ClS db result, 



6S 



(^'°)(/(x),/'(x))-^[^(°'^)(/(x),r(x))]. 



Ti{y) 



(A.40) 



(A.41) 



It is interesting to note that this is almost an expression of the Euler-Lagrange equation 
for a single dynamical variable. To see this, note that the derivative may more succinctly — 
though less explicitly — be written as 



5^_dg_ d_ 
6f df dx 



dg_ 
df 



(A.42) 



where by |j and ^ we mean the derivatives of g with respect to the first and second arguments 
respectively, evaluated at (/(x), /'(x)). Making the suggestive replacements g ^ L, x t, 
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Table A.l: Some useful functional derivatives. The notations and g^'^'^'> mean derivatives 
of g with respect to its first and second arguments respectively. 

Functional S Functional derivative -At-t 
^fm 

J dygUm g'UXx)) 

j dyg{y)f{y) -g\x)f{x) 

I dyg{f{y)J'{y)) g^^''\f{x)J'{x))-^[g^^'^\f{x)J'{x))] 

f{y) S{x-y) 



f ^ q, and q' = q, we have 



and setting 



6S 



dq 



d_ 

di 



dq 



(A.43) 



— yields the Euler-Lagrange equation in its usual form. 
Table A.l presents a summary of several useful functional derivatives, most arising as special 
cases of Eq. (A. 41). 



Appendix B 

Details of two-dimensional simulations 



B.l Simulation using the PGPE 



Here we outline our procedure for determining the properties of the C region and the steps 
used to create initial states for the PGPE solver. The C region itself is characterised by the 
cutoff momentum while the initial states are characterised by the energy Ec and number 
Nq. We want to obtain values of these three properties that are consistent with a specified 
temperature T and total number of atoms A^. 

B.1.1 Hartree-Fock-Bogoliubov analysis 

To generate an initial estimate of the C region parameters we solve the self-consistent Hartree- 
Fock-Bogoliubov (HFB) equations in the so-called Popov approximation^ [58] to find an ap- 
proximate thermal state for the system at a temperature T. The resulting state is a Bose 
Einstein distribution of quasiparticles interacting only via the mean-field, expressed in terms 
of the quasiparticle amplitudes Wk and Wk- 

Occupations for the C region field may be computed directly from the quasiparticle occu- 
pations via 



where Nb is the Bose-Einstein distribution and is the quasiparticle energy that is obtained 
by solving the Bogoliubov-de Gennes equations self-consistently [58]. This allows us to compute 
the cutoff as the maximum value of ||k|| consistent with sufficient modal occupation: 




We choose ricut = 5 for the sufficient occupation condition on the C region modes. 

^It is amusing to note the commentary in Rcf. [126] that maintains that Popov himself never suggested this 
trick, and would only have deemed it valid very near the transition temperature! 




(B.l) 
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The number of atoms below the cutoff may be computed directly from the sum of the 
condensate number A'^o and the number of C region excited state atoms, A'^ic: 

A^c = iVo + A^ic, where N^c = J] ^k- (B.3) 

kGC\{0} 

For the total energy below the cutoff, we use the expression 



k6C\{0} 



where A = 5'(A'o + 2A''ic). Rearranging, this is 



E^ = ^2 i^c + N!c) + E E^[NBm-vl]. (B.5) 

kec\{o} 

The expression in Eq. (B.5) differs from Eq. (22) of Ref. [58] as we have retained the zeroth 
order (constant) terms that are required to match the energy scale of the HFB analysis to the 
zero point of energy in the classical field simulations. 



B.1.2 Initial conditions for fixed total number 

A simple comparison between simulations at varying temperatures can only be carried out 
if the total number of atoms is fixed. This presents a problem in our simulations: although 
the number of atoms and energy of the C region can be directly specified (see section B.1.3), 
we may only determine the total number after performing a simulation. This is because the 
number of atoms in the I region depends on the temperature and chemical potential that are 
calculated by ergodic averaging of the C region simulations. 
Formally, this may be stated as a root finding problem: solve 

A^(A^c,^c) = iVtot (B.6) 

with initial guess provided by the solution to the HFB analysis in section B.1.1. Although both 
A'^c and Ec affect the total number A^, we choose to fix A'^c to the initial guess and to vary Ec 
until the desired total number is found. 

We note that evaluating the function N{Nc;, Ec) is very computationally expensive and 
difficult to fully automate because it involves a simulation and several steps of analysis. For 
this reason we use a nonstandard root finding procedure: For the first iteration we simulate 
three energies about the initial guess Eq such that the results crudely span Ntot] these three 
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simulations can be performed in parallel which significantly reduces the time to a solution. A 
second guess was obtained by quadratic fitting of Ec as a function of which gives accurate 
to within about 5% of Atot- An additional iteration using the same interpolation method takes 
A^ to within 0.3%, which we consider sufficient. 

We note that changing Eq during the root finding procedure means we have no direct 
control over the final temperature of each specific simulation. In our case this is not a problem 
because we only require a range of temperatures spanning the transition. In principle one could 
solve for a given temperature by allowing A'c to vary in addition to Ec- 



B.l. 3 Initial conditions for given Ec and A^c 

We compute initial conditions for the C region field in a similar way to Ref. [31]. Using the 
representation for the C region given by Eq. (3.9), the task is to choose appropriate values for 
the {cn}. As a first approximation, choose the smallest value for a momentum cutoff K' such 
that the field with coefficients 



has energy greater than Ec- Here A is chosen so that the field has normalisation corresponding 
to A'c atoms, and 6^ is a randomly chosen phase which is fixed for each mode at the start of 
the procedure. The random phases allow us to generate many unique random initial states at 
the same energy. 

By definition, the field defined by Eq. (B.7) has energy slightly above the desired energy. 
This problem is solved by mixing it with the lowest energy state: 




for < ||k|| < K', 
for |k| > K', 



(B.7) 




for n = 



elsewhere. 



(B.8) 



using a root finding procedure to converge on the desired energy Ec- The scheme generates 
random realisations of a non-equilibrium field with given Ec and A'c which are then simulated 
to equilibrium before using ergodic averaging for computing statistics. 
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B.2 I region integrals 



Our assumed self-consistent Wigner function (section 3.3.2) for the I region atoms takes a 
particularly simple form in the homogeneous case: 

The above-cutoff density may then be found by direct integration: 

ni{x)= [ d'^kWi{k,x), (B.IO) 



X2 



In 



(B.ll) 



In a similar way, the assumed Wigner function allows any desired physical quantity to be 
estimated via a suitable integral. A particular quantity of interest in the current work is the 
first-order correlation function, which can be obtained from the Wigner function as [91] 

gS'^ (x, x') = / d^k e-^^'^-^''-''') Wi (k, ^) . (B. 12) 



This integral is of the general form 

/i(r) = / d'k (B.13) 

^||k||>A' <2 i 

for constants A and B. Noting that Ji depends only on the length r of ||r||, and transforming 
k to polar coordinates (k, 9), we have 

oo , /•2tt 

K 



K e ■ - ± Jo 



h{v)=l dK-^-,^^1 rf^e-™^ (B.14) 



/oo 
^ ^ %^.^+B_l 2[r(|)]VoM, (B.15) 

(see Ref. [55, p902] for the Bessel function identity). 

Thus we obtain ^[^^(x, x') in terms of a one- dimensional integral which may be performed 
numerically: 

G(1)(x,x') = - r dK ^ • (B.16) 
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B.3 Vortex detection 

The defining feature of a "cliarge-m" vortex is tliat tlie phase 9 of the complex field ^/'(x) = 
|'?/'(x)|e*^*-'''' changes continuously from to 2m7r around any closed curve that circles the vortex 
core. We express our field ip on a, discrete grid in position space; the aim of vortex detection 
is then to determine which grid plaquettes (that is, sets of four adjacent grid points) contain 
vortex cores. 

To obtain the phase winding about a plaquette, first consider the phase at two neighbouring 
grid points A and B. We are interested in the unwrapped phase difference A^ab between the 
grid points; unwrapping ensures that the phase is continuous between A and B. (In the discrete 
setting such continuity is poorly defined; the best we can do is to correct for the possibility of 
2ti phase jumps by adding or subtracting factors of 2tx so that |A6'ab| < tt.) The unwrapped 
phase differences around a grid plaquette tell us a total phase change ^wrap = Tli A6'i,i+i = 2m7r 
where m G Z is the winding number or "topological charge" . 

Due to the necessity of unwrapping the phase, a four-point grid plaquette cannot unam- 
biguously support vortices with charge larger than one. Luckily, such vortices are energetically 
unfavourable in 2D Bose gases [100, page 83] so we need only concern ourselves with detecting 
vortices with winding number ±1 in this work. The positions obtained from a given run of our 
vortex detection algorithm are the labelled {r^} and {r^} for winding numbers +1 and — 1, 
respectively. 



Appendix C 

An imaginary time method with fixed 

normahsation 



Ground states for the GPE are often found using the so-called "imaginary time" method [21]. 
The prescription is to replace the time t with the "imaginary time" r = it and the potential 
V with y — yU in the equation of motion. One then selects a desired chemical potential /i, and 
evolves the new equations in r until a stationary state 0o is reached. This stationary state 
corresponds to a ground state of the GPE. 

A primary advantage of the imaginary time technique is the extreme simplicity of implemen- 
tation, given the existence of a numerical solver for the real-time version of the equation. An 
important disadvantage is that the normalisation is controlled only indirectly via the chemical 
potential, and it is common to desire a fixed total number of atoms in a simulation instead. In 
the following we describe a method of achieving this by evolving /i along with the imaginary 
time evolution of the fields. 



C.l The usual imaginary time method 

We start by investigating the reasons why the traditional imaginary time method works. In 
the case that the interaction strength is zero, we have the linear Schrodinger equation and the 
reasoning behind the method is straightforward: The coefficient of the lowest energy eigenstate 
evolves as e~^°'^/^^ and therefore becomes large compared to coefficients e~^^'^^^^ as r — )■ oo, 
where Ek is the energy of the kih eigenstate. The wavefunction converges exponentially to the 
ground state as a result, provided that it has nonzero overlap with the initial state. 

In the case that the Hamiltonian is nonlinear there is no simple formal solution to rely on, 
so the justification for the imaginary time method must be modified. The appropriate starting 
point is to realise that we seek to minimise the Hamiltonian under the constraint that the 
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number of atoms, is fixed. That is, we wisli to solve the constrained minimisation problem, 

minimise T-Oid] 

^ ^ (C.l) 
subject to U[(j)\ = N 

where A/'[0] = / dV \(f)f is the normalisation functional and N is the desired normalisation. 
Proceeding via Lagrange multipliers, we may define 

/C[0]=H[0]-/i(Ar[0]-iV), (C.2) 

and solve for the critical point of /C. It is important to note that the critical point of the 
Lagrange function is not a minimum with respect to the variables {0, /i}; it is a saddle point. 
To find the critical point we must solve the two equations 

The first of these is the condition that the right hand side of the GPE with an extra energy 
offset /i is zero, while the second recovers the number constraint: 

(iJsp-/i)0 + t/|0p0 = O and 7V[0]-A = O. (C.4) 

In the conventional imaginary time method, we solve a restricted version of this problem by 
fixing yU and minimising /C. This works because the critical point of JC with respect to is a 
minimum for any fixed /i, even though it is not a minimum when the combined variables {(f), fi} 
are considered. (To see this, note that JC is also a valid Hamiltonian with an energy offset //.) 

Minimisation of /C with fixed fi may be achieved by continuous-time steepest-descent optimi- 
sation and this corresponds to the imaginary time evolution. This may seem like a coincidence 
but the fact that imaginary time evolution minimises the energy is a generic property of Hamil- 
tonian systems in complex phase space. 

To see this, note that the GPE may be written 

On the other hand, the steepest descent direction for /C is —j^, so an evolution equation which 
takes us toward the ground state is 

The factor of h on the left is included only to emphasise the similarities with Eq. (C.5). Taken 
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together, this shows why the imaginary time method works: it corresponds to a steepest-descent 
minimisation of the Hamiltonian, with chemical potential added as a Lagrange multiplier for 
the number constraint. 



C.2 The imaginary time method with varying fi 

The conventional imaginary time method as described above does not allow us to satisfy the 
normalisation constraint directly; to do this we must modify fi during the evolution. Surpris- 
ingly, the appropriate evolution law turns out to be the uphill evolution 

^ - +— (C 7) 

Combined with Eq. (C.6) this takes us to the appropriate saddle point of K. representing the 
solution of Eqs. (C.3). 

A discrete-time version of this saddle point finding algorithm is described in [82, §14.2 p. 429] 
and the local convergence of the continuous version proven in Ref. [6, Theorem 1]. For con- 
vergence to hold, the Hessian of /C with respect to (p must be positive definite at the critical 
point. We note that this is the case for our system, because for every fixed fi the critical point 
is in fact a minimum with respect to variations of 0. While the convergence theorem holds 
only locally, we find that the Thomas-Fermi solution is a sufficiently accurate starting point for 
reliable convergence in practice. 

We finish by giving some intuitive reasoning behind the saddle point algorithm. To simplify 
the exposition, we work with an analogous minimisation problem over R"': 

minimise f(x) 
subject to g[x) = 

where f,g: R" — !■ R. The Lagrangian function for this minimisation is 

£(x,/i) = /(x)-/i(7(x), (C.9) 

and we wish to solve for a critical point (xq, /io) of C In nearly all cases (xq, /io) is a saddle 
point because to first order near the critical point, 



fxg[x) = /i 



■ (x - Xo) + . . . 



(C.IO) 



/iV(7(xo) ■ (x-xo), (C.ll) 
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where the constant term disappears because the constraint g is satisfied at the critical point. 
This first-order approximation therefore amounts to a sum of saddle-hke terms of the form 

— Xo,i). 

In the case that the principle directions of the saddle point correspond to the coordinate 
axes, for example 



one can intuitively expect that a steepest descent in x and steepest ascent in y would converge 
to the saddle point. In the Lagrangian case, the principle directions do not line up with the 
coordinate axes, and instead we have a saddle point intuitively similar to the function 




(C.12) 




(C.13) 



Nevertheless, the evolution 




(C.14) 



dC 



(C.15) 



still converges to the saddle point in many cases, provided we start sufficiently close. As noted 
above, this happens when the Hessian of C with respect to x is positive definite [6, Theorem 
1]. This is the case when finding ground states of the GPE using the saddle point algorithm. 



Appendix D 

Additional work 



The following two papers contain additional work completed during the course of the PhD, but 
which do not tie into the main theme of the thesis. 

In the first paper we report the first discovery of a Bell inequality for observables with 
continuous spectra. In this paper the current author contributed the no-go proof on pages 
three and four, in collaboration with E. G. Cavalcanti: There are no local hidden variable 
inequalities possible if one considers only the first-moment correlations between continuous 
variables at different sites. The current author also helped with the final stages of writing the 
paper. 

The second paper proposes a computationally efficient alternative to the traditional Uhlmann- 
Jozsa fidelity measure between two quantum states. The current author contributed section IV 
that discusses computational efficiency, including the optimised C implementations, and addi- 
tional optimisation of the slower Matlab codes. The current author also helped with polishing 
of the manuscript whole. 

Note for the arXiv version: Verbatim copies of the two papers, Refs. [18] and [84] were 
attached here in the version of this thesis submitted for examination. These have been omitted 
from the arXiv version due to technical constraints. 
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